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An equation of motion for relativistic compact binaries is derived through the third post- 
Newtonian (3PN) approximation of general relativity. The strong field point particle limit and 
multipole expansion of the stars are used to solve iteratively the harmonically relaxed Einstein 
equations. We take into account the Lorentz contraction on the multipole moments defined in our 
previous works. We then derive a 3PN acceleration of the binary orbital motion of the two spherical 
compact stars based on a surface integral approach which is a direct consequence of local energy 
momentum conservation. Our resulting equation of motion admits a conserved energy (neglecting 
the 2.5PN radiation reaction effect), is Lorentz invariant, and is unambiguous: there exist no un- 
^~>) ■ determined parameters reported in the previous works. We shall show that our 3PN equation of 

' motion agrees physically with the Blanchet and Faye 3PN equation of motion if A = —f 987/3080, 

, where A is the parameter which is undetermined within their framework. This value of A is consis- 

_ ' tent with the result of Damour, Jaranowski, and Schafer, who first completed a 3PN iteration of 

' the ADM Hamiltonian in the ADMTT gauge using dimensional regularization. 
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I. INTRODUCTION 
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0^ ' Tremendous effort has been made for direct detection of gravitational waves. Interferometric gravitational wave 
detectors such as GEO600 [1], the Laser Interferometric Gravitational Wave Observatory(LIGO) [2], and TAMA300 
[3] have been operating successfully. They have been actively investigating the data and reported results [4-9]. 
One promising source of gravitational waves for those detectors is a relativistic compact binary system in an 
' inspiraling phase. The detectability and quality of measurements of astrophysical information of such gravitational 
wave sources rely on the accuracy of our theoretical knowledge about the waveforms. A high-order, say, third- or 
fourth-order, post-Newtonian equation of motion for an inspiraling compact binary is one of the necessary ingredients 
to construct and study such waveforms [10-13]. 
I In addition to the purpose of making accurate enough templates, the high-order post-Newtonian approximation 

^jPj for an inspiraling compact binary is a fruitful tool, for example, to construct astrophysically realistic initial data for 
• • numerical simulations [14-17] and to estimate innermost circular orbits [18,19]. 
. ^ ' The post- Newtonian (PN) equation of motion for relativistic compact binaries, a fundamental tool employed for 
, the above purposes, has been derived by various authors (see reviews, e.g., [20-22]). The equation of motion for 
\^ ' a two point-masses binary in a harmonic coordinate up to 2.5PN order, at which the radiation reaction effect first 
appears, was derived by Damour and Deruelle [23,24] based on the post-Minkowskian approach [25]. These works used 
Dirac delta distributions to express the point-masses mathematically, therefore they inevitably resorted to a purely 
mathematical regularization procedure to deal with divergences arising due to the nonlinearity of general relativity. 
Damour [20] gave a physical argument known as the "dominant Schwarzschild condition" which supports the use of 
Dirac delta distributions and their regularization up to 2.5PN order. 

Direct validations of the Damour and Deruelle 2.5PN equation of motion have been obtained by several works 
[26-30]. Grishchuk and Kopeikin [26] and Kopeikin [27] worked on extended bodies with weak internal gravity. 
On the other hand, the present author, Futamase, and Asada derived the 2.5PN equation of motion [29] using the 
surface integral approach [31] and the strong field point-particle limit [32]. All the works quoted above agree with 
each other. Namely, our work [29] shows the applicability of the Damour and Deruelle 2.5PN equation of motion 
to a relativistic compact binary consisting of regular stars with strong internal gravity. We mention here that stars 
consisting of relativistic compact binaries, for which we are searching as gravitational wave sources, have strong 
internal gravitational field, and that it is a nontrivial question whether a star follows the same orbit regardless of the 
strength of its internal gravity. 
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A 3PN iteration was first reported by Jaranowski and Schafer [33]. There a 3PN Arnowitt-Deser-Misner (ADM) 
Haniiltonian in the ADM Transverse Traceless (ADMTT) gauge for two point-masses expressed as two Dirac delta 
distributions was derived based on the ADM canonical approach [34,35]. However, it was found in [33] that the 
regularizations they had used caused one coefficient Wkinetic undetermined in their framework. Moreover, they later 
found another undetermined coefficient in their Hamiltonian, called ecstatic [36]. Origins of these two coefficients 
were attributed to some unsatisfactory features of regularizations they had used, such as violation of Leibniz's rule. 
The former coefficient was then fixed as Wkinetic = 41/24 by a posteriori imposing Poincare invariance on their 3PN 
Hamiltonian [37]. As for the latter coefficient, Damour et al. [38] succeeded in fixing it as ^static = 0, adopting 
dimensional rcgularization. Moreover, with this method they found the same value of Wkinetic as in [37], which ensures 
Lorentz invariance of their Hamiltonian. 

On the other hand, Blanchet and Faye have succeeded in deriving a 3PN equation of motion for two point-masses 
expressed as two Dirac delta distributions in a harmonic coordinate [39,40] based on their previous work [28]. In 
this approach, they have assumed that two point-masses follow regularized geodesic equations (more precisely, they 
have assumed that the dynamics of two point-masses is described by a regularized action, from which the regularized 
geodesic equations were shown to be derived). The divergences due to their use of Dirac delta distributions were 
systematically regularized with the help of generalized Hadamard's partie finie rcgularization, which they have devised 
[41] and carefully elaborated [42] so that it respects Lorentz invariance. They found, however, that there exists one 
and only one undetermined coefficient (which they call A). (But see [43] for the recent achievement of Blanchet and 
his collaborators.) 

Interestingly, the two groups independently constructed a transformation between the two gauges and found that 
these two results coincide with each other when there exists a relationship [44,45] 



Therefore, it is possible to fix the A parameter from the result of [38]. However, applicability of mathematical 
regularizations to the current problem is not a trivial issue, but an assumption to be verified, or at least supported 
by convincing arguments. There is no argument such as the "dominant Schwarzschild condition" at the 3PN order. 
Thus, it seems crucial to achieve unambiguous 3PN iteration without introducing singular sources and to support (or 
give counterevidence against) the result of [38]. 

In this paper, we derive a third post-Newtonian equation of motion in a harmonic coordinate applicable to inspiraling 
binaries consisting of regular spherically symmetric compact stars with strong internal gravity, using yet another 
method which is based on our previous papers ( [46] and [29], referred to as paper I and paper II, respectively, 
henceforth). Namely, to treat strong internal gravity of the stars, we have used the strong field point-particle limit 
and surface integral approach. This point-particle limit enables us to incorporate a notion of a self-gravitating point 
particle into general relativity without introducing singular sources. 

In [47], we made a short report on our results. In this paper, we shall present the full explanation of our method 
and show our results. 

The organization of this paper is as follows. In Sees. II and HI, we briefly explain the basics of our method (see 
papers I and II for more details). After a short explanation on the structure of the 3PN equation of motion, we then 
compute the 3PN gravitational field required to derive the 3PN mass-energy relation from Sees. IV to VI. The 3PN 
mass-energy relation and the 3PN momentum-velocity relation are then derived in Sees. VII and VIII. Section IX 
describes how we evaluate the 3PN gravitational field necessary to derive a 3PN equation of motion. In these sections, 
we write down the formulas we have used and show several examples of our computations only, since the intermediate 
results are generally too lengthy to write down. We then show a 3PN equation of motion we obtained in Sees. X and 
XI and we shall compare it with the Blanchet and Faye 3PN equation of motion in Sec. XII 1. There, we found that 
our result is physically consistent with their result with A — —1987/3080. Section XII 2 is devoted to a summary of 
our formalism and discussion. Some useful formulas and supplementary computations are shown in Appendixes. 

Throughout, we use units where c = 1 = G. A tensor having alphabetical indexes, such as or x, denotes a 
Euclidean three-vector. We raise or lower its indexes with a Kronecker delta. An object having Greek indexes, such 
as a;**, is a four- vector. Its indexes are raised or lowered by a flat Minkowskian metric. 



In this section, we briefly review the strong field point-particle limit and associated scalings of the matter on the 

initial hypersurface. See [29,32,46,48 52] for more details. 

We first introduce a nondimensional small parameter e which represents slowness of a star's typical orbital velocity 
^orb ^"^^ thMS which is a post-Newtonian expansion parameter. 
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II. STRONG FIELD POINT-PARTICLE LIMIT AND SCALINGS 
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where we set = dx'' /dr of order unity. We call r the Newtonian dynamical time [32]. t is a natural dynamical time 
scale of the orbital motion [32,49]. 

Then the post-Newtonian scaling implies that the typical scale of the mass of the star scales as e^. 

Henceforth, we call the {t,x^) coordinate the near zone coordinate. 

A. Strong field point-particle limit 

One would think that a point-particle limit may be achieved by setting the radius of the star to zero. In general 
relativity, however, by this procedure the star cannot become a "point-particle" , rather it becomes an extended body 
(black hole) whose radius is of order of its gravitational radius. One solution to avoid this conceptual problem was 
proposed by Futamase [32]. Following Futamase, we achieve a point-particle limit by letting the radius of the star 
shrink at the same rate as the mass of the star. This limit is nicely fit into the post-Newtonian approximation, since 
from the post- Newtonian scaling of the mass (0(e2)), the radius of the star is ©(e^) ^ as we take the post-Newtonian 
limit (e ^ 0). This point-particle has finite internal gravity since a typical scale of the self-gravitational field of the 
star, the mass over the radius, is finite. Thus we call this limit the strong field point-particle limit. 

B. Surface integral approach and body zone 

We construct a sphere for each star, called the body zone Ba for the star A {A = 1, 2) [32], which surrounds each 
star and does not overlap the other. More specifically, the scalings of the radius and the mass of the star motivate us 
to introduce the body zone of the star A, Ba = {x'^\\x — za{t)\ < ei?A} and a body zone coordinate of the star A, 
a\ = e~^(a;* — Za{t)). Here Za{t) is a representative point of the star A, e.g., the center of mass of the star A. Ra, 
called the body zone radius, is an arbitrary length scale (much smaller than the orbital separation while cRa is larger 
than the radius of the star for any e) and constant (i.e., dRA/dr = 0). Using the body zone coordinate, the star does 
not shrink when e ^ 0, while the boundary of the body zone goes to infinity. Thus, it is appropriate to define the 
star's characteristic quantities such as a mass using the body zone coordinate. 

On the other hand, the body zone serves us as a surface OB a, through which gravitational energy momentum fiux 
flows and in turn it amounts to gravitational force exerting on the star A. Since the body zone boundary OBa is 
far away from the surface of the star A, wc can evaluate explicitly the gravitational energy momentum flux on OBa 
with the post-Newtonian gravitational fleld. After computing the surface integrals, we make the body zone shrink to 
derive the equation of motion for the compact star. 

Effects of the internal structures of the compact stars may be coded in, e.g., multipole moments of the stars. These 
moments in turn appear in the gravitational energy momentum flux in the surface integral approach and affect the 
orbital motion. 



C. Scalings on initial hypersurface 

Following the works [50-52], we take the initial value formulation to solve Einstein equations. As the initial data 
for matter variables and gravitational field, we take a set of nearly stationary solutions of the exact Einstein equations 
representing two widely separated fluid balls. We assume that these solutions are parametrized by e and the matter 
and the field variables have the following scalings on the initial hypersurface. 

The matter density scales as e^^ (in the (t,x^) coordinate), implied by the scalings of the mass and the radius of 
the star. The internal time scale of the star is assumed to be comparable to that of the binary orbital motion. We 
assume that the star is pressure-supported. 

Prom these initial data we have the following scalings of the star ^'s stress energy tensor components in the body 

zone coordinate, T^^: = 0(e~^), = 0{e~'^), — 0{e^^) on the initial hypersurface. Here the underlined 
indexes mean that for any tensor S", S- = e~^5*. In paper I, we regard the star A's body zone coordinate as a 
Fermi normal coordinate of the star and we have transformed Ty" , the components of the stress energy tensor of the 
matter in the near zone coordinate, to T^" using a transformation from the near zone coordinate to the Fermi normal 
coordinate at IPN order. It is diflacult, however, to construct the Fermi normal coordinate at a high post-Newtonian 
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order. Therefore, we shall not use it. We simply assume that for T'^" (or rather for A^'', the source^ term of the 

harmonically relaxed Einstein equations), TJ/ = 0(e~^),T^- = 0(e~^), and = 0{e~^), as their leading scalings. 

As for the field variables on the initial hypersurface, we simply make a reasonable assumption that the field is of 
2.5PN order except for the field determined by the constraint equations. If we take random initial data for the field 
[50] supposed to be of IPN order, they are irrelevant to the dynamics of the binary system up to the radiation reaction 
order [52]. Thus, we expect that the 2.5PN order free data of the gravitational field on the initial hypersurface do not 
affect the orbital motion up to 3PN order. 



III. FORMULATION 
A. Field equation 

As discussed in the previous section, we shall express our equation of motion in terms of surface integrals over 
the body zone boundary where it is assumed that the metric slightly deviates from the flat (auxiliary) metric 7]^^" = 
diag(— e^, 1, 1, 1) (in the near zone coordinate (r, x')). We define a deviation field /i'*'' as 

/i^- = - (3.1) 

where g is the determinant of the metric. Indexes of h^" are raised or lowered by the flat metric. 

Now we impose a harmonic coordinate condition on the metric ft-^^.^ = 0, where the comma denotes a partial 
derivative. In the harmonic gauge, we can recast Einstein equations into a relaxed form, 

□ /i"^ = -167^A'^^ (3.2) 

where □ = 'q'^^dfj^di, is the flat d'Alembertian and 

Ki"" = + X^"°''^,oci3, (3.3) 
e^- = (_5)(T''- + if^), (3.4) 

X^""''^ = -^i'^"''^'^^ - h°''^h^''). (3.5) 

Here, T^^"^ and f^"^ denote the stress-energy tensor of the stars and the Landau-Lifshitz pseudotensor [53]. In consis- 
tency with the harmonic condition, a local energy momentum conservation law is expressed as 

A^^t. = 0. (3.6) 

Note that x^""'^ ,a(3 itself is divergence- free. 

Now we rewrite the relaxed Einstein equations into an integral form. 



,3Ai'"{T-c\x-!j\.!/'^: 



ft^^(T,x^) = 4/ d'y ^ ,J-- ' ' +h^jj'{T,x^), (3.7) 

Jc{t,x>') f - y\ 

where C(r, a;^) means the past light cone emanating from the event (r, x'^). C{t, x'') is truncated on the r = initial 
hypersurface. h'^ is a homogeneous solution of a homogeneous wave equation in flat spacctime. In this paper, we 
shall adopt the no-incoming radiation condition (see, e.g., [54]) by taking sufficiently large r, i.e., h'^ = 0. 

We solve the Einstein equations as follows. First we split the integration region into two zones: the near zone and 
the far zone. 

The near zone is the neighborhood of the gravitational wave source where the wave character of the gravitational 
radiation is not manifest. In this paper, we define the near zone as a sphere centered at some fixed point, enclosing 
both of the stars, and having a radius TZ/e, where K is arbitrary but larger than the binary separation and the 
gravitational wavelength. The scaling of the near zone radius is derived from the e dependence of the wavelength of 
the gravitational radiation emitted by the binary due to its orbital motion. The center of the near zone sphere would 
be determined, if necessary, for example, to be the center of mass of the near zone. The outside of the near zone is 
the far zone where the retardation eff'ect of the fleld is crucial. 

We evaluated the integrals over the far zone which contribute to the near zone field where the stars reside using the 
direct integration of the relaxed Einstein equations (DIRE) method, which was initiated by Will and his collaborators 
[30,55,56]. DIRE directly and nicely fits into our formalism since it utilizes the relaxed Einstein equations in the 
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harmonic gauge. Although wc do not show our cxphcit computation in this paper, we have followed the DIRE 
method and checked that the far zone contribution does not affect the equation of motion through 3PN order. In 
fact, Blanchet and Damour [57], and later Pati and Will [56], showed that the far zone contribution to the near zone 
field affects (physically) the orbital motion starting at 4PN order. Henceforth we shall focus our attention on the near 
zone contribution to the near zone field, 

/.-(r,.^) = 4 / d3,^^£(L^4^Z^.^, (3.8) 

Jn \x-y\ 

where denotes the near zone and we attach the subscript N to h.^" to clarify that they are quantities in the near 
zone. 



B. Near zone contribution 



For the near zone contribution, we first make retardation expansion and change the integral region to a r = const 
spatial hypersurface 

,,,^=^Y^tl}lfl_Y f d'y\g-y\n-^A%^{r,y'';e). (3.9) 

Note that the above integral depends on the arbitrary length TZ in general. The cancellation between the 7?.-dependent 

terms in the far zone contribution and those in the near zone contribution was shown by Pati and Will [56] throiigh 
sufficient post-Newtonian order. Moreover, the findings in [56,57] make us expect that we can remove 7?,-dependent 
terms via a suitable gauge transformation and then take a formal limit of infinite TZ in an equation of motion up to 
the 3PN level. With this expectation in mind, in the following we shall omit the terms with negative powers of TZ 
{JZ~^; k > 0) which vanish in the TZ infinite limit, while we shall retain terms with positive powers of TZ {TZ^\ k > 0) 
and logarithmic terms (In 7?.) and see as a good computational check that these 7?.-dependent terms can be gauged 
away from our final result. Another reason to keep logarithmic terms {InTZ) is to make the arguments of all possible 
logarithmic terms nondimensional. 

Second wc split the integral into two parts: contribution from the body zone B = Bi U B2, and from elsewhere, 
N/B. Schematically we evaluate the following two types of integrals (we omit the indexes): 

h = Y,{hBn + hN/Bn), (3-10) 

>.s, = 4tlll J2 ! <P„./Mi±f!fe2, ,3.U) 



n\ \dTj \rA-e'^dA\^ " 

^ ^ __,_d y , (3.12) 



n\ \dTj J 1^-^^-"' 

where fA = x — ZA and /(r, x'') is some function. We shall deal with these two contributions successively in the 
following. 



1. Body zone contribution 



As for the body zone contribution, we make a multipole expansion being concerned with the scaling of the integrand, 
i.e., A'^'' in the body zone. For example, the n = terms in Eq. (3.9), h'^^^^, give 
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/ T3T r\k „k OT<kt> k „l tiT<klm>„k„l ~,m\ 

+0(ei^), (3.13) 

/pi jki^k Qj<kl>i k„l \ 

/^F„=o=4.^ E (^ + ^^^+^^ (3-14) 



A=1.2 




ykij k oy<kl>ij k„l r.y<klm>ij k „l „7yi 

2^A 'A , A^'^yl 'a'A , fi^^A 'A'A'A 



2ri 2r^ 



(3.15) 



where ta = |'^4|- The quantity with <> denotes symmetric and trace-free (STF) on the indexes between the brackets. 
To derive the 3PN equation of motion, we need /i"^^ up to 0(e^°) and /i^' and /i'^ up to O(e^). 
In the above equations we defined multipole moments of the star A as 



ij^ —e 



[ d^aAA^Ja^, (3.16) 

Jba 

Jf ^ = e^ [ SaAlC^c^, (3.17) 

JBa 

Zf'^^e'f dW|af, (3.18) 
Jba 

where the capital index denotes a set of collective indexes, 7; = 11^2 ■ ■ ■ ii and cr^ = ar^a'J ■ ■ ■ a\. Then = I^, 

= IJ^ , and = . Wc simply call the four-momentum of the star A, P\ the (three-)momentum, and 
the energy. Also we call D\ the dipole moment and the quadrupole moment. 

Then we transform these moments into more convenient forms using the conservation law Eq. (3.6). In the following, 
ovcrdot denotes r time derivative, and yA = y — za- Noticing that the body zone remains unchanged (in 
the near zone coordinate), i.e., Ra = 0, we have 

P\ = Plv\ + Q\ + e'^, (3.19) 

Jl = \ [mJ + j + v^J^D^} + y^Qt, (3.20) 

+ e'Q<^v^2+e'Rf^ + le'^, (3.21) 
3 



k^3 - A^^'^\ (3.22) 



2 



where 



Ma = 2e* / d^aAa^A^, (3.23) 

Jba 



IBa 

= e-* / dS^ (A]v™ - v^A^^) yl^y\, (3.24) 
Joba 

Rf'^ ^ e-' / dSm (a7 - <Ay) yfyi,, (3.25) 

JdBA ^ ' 



and 



.l^^^' = ^Jfv^ + e^i;^^^'^) + + e^^:^. (3.26) 
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[ ] and ( ) denote antisymmetrization and symmetrization on the indexes between the brackets. M^J is the spin of the 
star A and Eq. (3.19) gives a momentum- velocity relation. Thus our momentum- velocity relation is a direct analogue 
of the Newtonian momentum- velocity relation. In general, we have^ 



jKii _ -j(Kii) 



21 

T+1 



J 



7Kiij 



+ 



21 

ITT" 



AK,-i[k,)i]j 



21 

T+T' 



7{Ki-i[k,)j]i 



(3.27) 
(3.28) 



and 



Kii 



2dlA 



l + l 



dr 



+ v\I 



(i tKi) 



Z + 1 



e Ha ' 



l + l 



dr 



l + l 



(3.29) 
(3.30) 



The surface integrals and i?^'*-' will be computed in Appendix A and do contribute to the field and the equation 
of motion starting at 3PN order. 



2. N/B contribution 

About the A''/B contribution, since the integrand A^^'' = —gt'^i^+x^'^"^^ ,a0 is at least quadratic in the small deviation 
field h^'^ , we make the post-Newtonian expansion in the integrand. Then, basically, with the help of (super)potentials 
g{x) which satisfy Ag{x) = f{x), A denoting a Laplacian, we have for each integral (e.g., the n = term in Eq. 
(3.12)) 



J N/B \x-y\ 



-47r5(x) 



dSk 



d{N/B) 



dg{y) 



1^ — ^ dy'' 



(3.31) 



We show a derivation of Eq. (3.31) in Appendix B. For n > 1 terms in Eq. (3.12), we use (super)potentials many 

times to convert all the volume integrals into surface integrals and the bulk terms ("— 47rg(.x)"). ^ 

Finding the superpotentials is one of the most formidable tasks, especially when we proceed to a high post-Newtonian 
order. Fortunately, up to 2.5PN order, all the required superpotentials are available [28,33]. At 3PN order, there 
appear many integrands for which we could not find the required superpotentials. To obtain a 3PN equation of 
motion, we use an alternative method. The details of the method will be explained later. 
Finally, we note that h^" = O(e^) as shown in paper II. 



C. General form of the equation of motion 

From the definition of the four-momentum, 

P^(r) ^e' f d^aAA^^^, 

JBa 

and the conservation law, Eq. (3.6), we have an evolution equation for the four-momentum, 

^ = / dSkA%^ + e-\\ I dSkA^T- 



(3.32) 



(3.33) 



^The equation in paper II corresponding to Eq. (3.30) has a misprint, though it does not affect the 2.5PN equation of motion. 
^Notice that when solving a Poisson equation Ag(x) = f{x), a particular solution suffices for our purpose. By virtue of the 
surface integral term in Eq. (3.31), it is not necessary to be concerned about a homogeneous solution of the Poisson equation. 
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Here we used the fact that the size and the shape of the body zone are defined to be time-independent (in the near 
zone coordinate), while the center of the body zone moves at the velocity of the star's representative point. 

Substituting the momentum- velocity relation, Eq. (3.19), into the spatial components of Eq. (3.33), we obtain the 
general form of the equation of motion for the star A, 

-4 



dBA JOB A 



+€- 



\JdBA JdBA J 



_dQ^_ fD^ 

dr ' dr2 ■ ^^-"^^^ 

All the right-hand side terms in Eq. (3.34) except for the dipole moment are expressed as surface integrals. We 
can specify the value of D\ freely to determine the representative point z\{t) of the star A. For example, we may 
call z\{t) corresponding to D\ = the center of mass of the star A from an analogy of the Newtonian dynamics. 

In Eq. (3.34), rather than the mass of the star A appears. Hence we have to derive a relation between the mass 
and PJ. We shall derive the relation by solving the temporal component of the evolution equation (3.33) functionally. 
In fact, at lowest order, we have shown in paper II that 

^ = 0(e2). (3.35) 

dr 

Then we define the mass of the star A as the integrating constant of this equation, 

mA = limPI. (3.36) 

e->0 

niA is the ADM mass that the star A would have if it were isolated. We took the e zero limit in Eq. (3.36) to ensure 
that the mass defined above does not include the effect of the companion star and the orbital motion of the star itself. 
Some subtleties about this definition are discussed in paper II. By definition, is constant. The procedure that we 
use to solve the evolution equation of and obtain the mass energy relation is achieved up to 3PN order successfully 
and the result will be shown later. 

Since the general form of the equation of motion is expressed in terms of surface integrals over the body zone 
boundary, we can derive an equation of motion for a strongly self-gravitating star using the post-Newtonian approx- 
imation. Effects of the star's internal structure on the orbital motion such as tidally induced multipole moments 
appear through the field and hence the integrand A^'' of the surface integrals. 



D. Lorentz contraction and multipole moments 



In this paper, we are concerned with a binary consisting of two spherically symmetric compact stars. In other 
words, all the multipole moments of the star defined in an appropriate reference coordinate where effects of its orbital 
motion and the companion star are removed (modulo, namely, the tidal effect) vanish. We adopt the generalized 
Fermi coordinate (GFC) [59] as a star's reference coordinate for this purpose 

Then a question specific to our formalism is whether the differences between the multipole moments defined in Eqs. 
(3.16), (3.17), and (3.18) and the multipole moments in GFC give purely monopole terms. This problem is addressed 
in Appendix C and the differences are mainly attributed to the shape of the body zone. The body zone Ba which is 
spherical in the near zone coordinate (NZC) is not spherical in the GFC mainly because of a kinematic effect (Lorentz 
contraction). To derive a 3PN equation of motion, it is sufficient to compute the difference in the STF quadrupole 
moment up to IPN order. The result is 

XT<i3> = T<iJ> _ r<y> _ 2 '^'^A <i j> , ^/ 3\ (o 
"^A — -'a.NZC -'A.GFC — ^^^^A^A + A {6.6() 

where Ia^^zc ^ -^A- ^agfc ^® ^^'■^ quadrupole moments defined in the generalized Fermi coordinate. 

As is obvious from Eq. (3.37), this difference appears even if the companion star does not exist. We note that we 
could derive the 3PN metric for an isolated star A moving at a constant velocity using our method explained in this 
section by simply letting the mass of the companion star be zero. Actually, SI^^-'^ above is a necessary term which 
makes the so-obtained 3PN metric the same as the Schwarzschild metric boosted at the constant velocity va in the 
harmonic coordinate. 
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E. On the arbitrciry constant Ra 



Our final remark in tiiis section is on the two arbitrary constants Ra- Since we introduce the body zones by hand, 
the arbitrary body zone radii Ra seem to appear in the metric, the miiltipole moments of the stars, and the equation 
of motion. As was discussed in detail in paper II, we proved that the surface integrals in Eq. (3.34) that we evaluate 
to derive the equation of motion do not depend on Ra through any order of the post-Newtonian iteration. As for the 
field and the multipole moments appearing in the field (mass, spin, quadrupoles, etc.), we reasonably expect that the 
i?^-dependent terms in the body zone contribution h'^ and N/B contribution h'^^^ cancel each other out, since the 
total field h^^"^ = h'^ + h^^g is obviously independent of Ra- In the cancellation, the multipole moments should be 
"renormalized" so that those moments do not depend on Ra ■ Such cancellation among ii^-dependent terms and the 
"renormalization" was demonstrated explicitly in paper I up to IPN order. 

Practically, the above observation enables us to discard safely all the ei?A-dependent terms except for logarithms 
of eRA- We retain In ei?yi-dependent terms to nondimensionalize the arguments of the logarithms. Thus, instead of 
performing renormalization of the multipole moments, we simply discard the eRA dependences other than the Inei?^ 
dependences. We emphasize here; tiiat we discard the ei?A-dependent terms in the field first and then evaluate the 
surface integrals in the general form of the equation of motion using the field which is independent of cRa- We then 
discard the ei?A-dependent terms arising in the computation of the surface integrals. 

Appendix D is devoted to an explanation on the renormalization of the stars' multipole moments. There, we also 
give a justification for our procedure of discarding through our computation of the field all the eii^-dependent terms 
other than the In eii^-dependent terms. 

IV. STRUCTURE OF THE 3PN EQUATION OF MOTION 

In the following sections, we shall derive an acceleration for two spherical compact stars through third post- 
Newtonian accuracy. 

First, we split the four-momentum, the dipole moment, and the Q\ integral into two parts, namely the part and 
the X part. 

P^e = I d^(^Ae%\ (4.1) 



Ba 



p/i — ^2 



e- / d^aAxT^,a?, (4.2) 

Jba 

D\q = e^ [ d^aAC^A^N, (4-3) 

-/Ba 

^e^ f d^aAa\xT^,a0, (4.4) 

JBa 

^ / dSk (e]v'= - <e]7) fA, (4.5) 

JOB A 

i dSu (xT^,<.p - v\xT^,<.p) V\. (4.6) 

JdBA ^ ' 



JdBA 

Correspondingly, we split the momentum- velocity relation Eq. (3.19) and the evolution equation for the four- 
momentum Eq. (3.33) into the 8 part and the x part. 

Now, and as well as can be expressed completely as surface integrals, and can be computed explicitly 
into functions of m^, va, and fi2. It is straightforward to compute them up to 3PN order, since we only need the 
2PN field to perform the surface integrals. The results are shown in Appendix E. There, we also compute the x 
part of the momentum- velocity relation and the evolution equation for Pa^- Then comparing these equations with 
Pax^ ^Ax' ^^'^ Qax^ found that the x part of these equations is an identity up to 3PN order. This observation 
then means that the nontrivial momentum- velocity relation and the mass-energy relation, an equation of motion, are 
obtained from the 6 part of Eqs. (3.19) and (3.33). 

Therefore, the equations that we have to evaluate to derive an evolution equation for the energy and an equation 
of motion are actually 
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\ dr y <3PN 

\ UT / <3pN 
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dr 



<2.5PN 



dSkwQ'^N + v'l f dSkloQ 

JdB, 



(4.7) 



= mi 



dvl 
d7 



<2.5PN 



dBi JdBi 



6 / dP^Q 
dr 

dr 



3PN 



vl+emrm-P^Q) 



r ,dv\ 
dr 



3PN 



dT2 ' 



(4.8) 



where for an equation or a quantity /, (/)<nPN and (/)nPN denote / up to riPN order and / at nPN order, respectively. 
<„/ and „/, on the other hand, denote an equation or a quantity / up to 0(e") and at 0(e"). In paper II, we found 
Q\q = 0{e^). For later convenience, in Eq. (4.8) we retain of order e"^, which appears at a 3PN equation 

of motion (in other words, we set <3D\q = 0). It should be understood that in the second line of Eq. (4.8), the 
acceleration dv\/dT should be replaced by the acceleration of an appropriate order lower than 2.5PN. We note that 
the X parts of Q^'* and i?^'*"* integrals and the multipole moments including Pj^^ and affect an equation of 



motion through the field and hence the integrands of Eqs. (4.7) and (4.8 
form of the 3PN equation of motion. 



Henceforth, we call Eq. (4.8) the general 



The explicit forms of the integrands loB^ = io[(— 5)*^!,] (o^^ dBA) are (see Appendix F), 



UTT,k 



+ ••• 
+ ••• , 



1677106*^ = \{5\5h + 5\5h - S^'Ski) {ih^^^'i 
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''+48/l"',.) 



84/1^ 



(4.9) 
(4.10) 

(4.11) 



The fields up to 2.5PN order, Kgh^'^, Krh^^, and <7h^^ , are listed in paper II. Thus, to derive the 3PN mass-energy 
relation and the 3PN momentum-velocity relation, we have to derive sh^^- To derive the 3PN equation of motion, we 
further need loh^'^ + ^h'^k- 

Up to 2.5PN order, the superpotentials required to compute the field could be found [28,33]. However, at 3PN 
order, it is quite difficult to complete the required superpotentials. We take another method to overcome this problem. 
We shall explain our method in detail in Sees. V and IX. Below, we begin our calculation by deriving 8/i^' and check 
(a part of) the 3PN harmonic condition <s,h^^ = 0. 



V. sW': N/B INTEGRALS 



In this section, we derive g/i^ 



<8' 



<8/iR + 4e' 



N/B 



\x — y\ dr"^ 



/ d^y\x-y\QK]i 

J N/B 



(5.1) 



where <8ft-B is the body zone contribution up to O(e^) and shown in Appendix A as Eq. (A9). The Q^'* integral 
contained in k^K^b found in Appendix A to vanish. We are concerned with the equation of motion for two spherical 
compact stars and we shall only retain monopole terms in <sh^. 

The second time derivative term in the retardation expansion (the last term in Eq. (5.1)) can be integrated explicitly 
via super-superpotentials (i.e., a particular solution of the Poisson equation with a superpotential as a source). The 
result is 



/ - 

Jn/b 



+ 



12 



6 \TZ/e 

+8P^p^ (^6\A,2 - d^.d,. + / 
np[Pi, /27^ 



+ 



(InS) 



■In 



+ (1^2), 



(5.2) 
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where /S.aa' = S'^-'d^j d j and d^i = d/dz\. The symbol (1 ^ 2) denotes the same terms but with the star's label 

' ^ A' ^ 

1 exchanged for 2. /C'^'^) is a superpotential satisfying A/^'"*^) = In 5, where S = ri + r2 + ri2, and its explicit 
expression is given in [33] as 

/('"^^ = 3^(-r-? + 3riri2 + rj^ - Snrj + Snsrj - r^) + ^(r^ - r^^ + 4) InS. (5.3) 

Now let us devote ourselves to an evaluation of the Poisson integral with the integrand s^]^- For this integrand, 

gA]^, it is difficult to find all the required superpotentials. We proceed as follows. First, wc split the integrand gA]^ 
into two groups, one whose members depend on the negative power or logarithms of S or both, and the other whose 
members do not. 

The S-dependent integrands are 

S - dependent parts of {-ih'''-^\h7 k,i + 2ih'^^^^\h'' k,i + 4/i*'=,r4/i^^,fe + 24/i^'',fe6/i^'*''''}- 

For the remaining integrands in gA]^, we found all the required superpotentials except for essentially two Poisson 
equations whose source terms are l/r\/r2 and (lnri)/r2. Then in the following, we split the integrands into three 

groups: 

(a) Direct-integration part = S-dependent group: 

the integrands which depend on inverse powers of S (5*^, k : negative integer), or have logarithms of S (InS), or 
both. 

(b) Superpotential part: 

the integrands which do not depend on inverse powers of S, nor have logarithms of S, and for which the required 
superpotentials are available, 

(c) Superpotential-in-series part: 

the integrands which do not depend on inverse powers of S, nor have logarithms of S, and for which the required 

superpotentials arc unavailable. 

We mention here that splitting the integrands into the S-dependent and the S-independent group is a rather rough 
technique; there may be some terms in the S-dependent group for which we could find superpotentials. We have made 
such a classification since for the S-dependent group it seemed difficult to find superpotentials. 

For example, let us take —4h^^''^ih'^k.i- Then 

[- Ah'^'''\h'^ kj]-QiY, = direct - integration part of {-Ah'^'''\h'^k,l} 
d,m\m2 j / 2 2 1 ^"12 ^2 



rlS \ riri2 rir2 ri2r2 rjr2 r^ri2 
16mlm2 , ( I ,^ ^ . ( I 2 2 ri2 2r: 



c2 '"i I {ni-vi)+\-— 2 2 '""2 1--2-^ (^i2-'yi) 

( r2 1 1 \ ^ , 
+ — 2 ('^2 • Vl) 

16mfm2 if 1 1 2 ^ _, ^ ,^ ■^iw-.-n 

"I To "2 ~~3 + ~3~ H 2~ + ~2 2~2 3 2~3~ i'^l ' ^l) 



2r2 




rl 


rlrl2 


r\r\2 


rM2 


2r2 


2rl 

4- 


rl 


nn2 


nri2 


r\r\2 



16m^m2 i f 2 1 ri2 r2 ~. ^ , , . , 

+ 52 rf r\T2 rlT2 nrf^ -2^2 +.3.2 +.3.3 ) • ^;i) 

+(1 ^ 2), (5.4) 



[-4/i'''''4/i'^fe,;] sp = superpotential part of {-^t^^''' k,i\ 

= {v\ - 3(ni • VY)n\) j-^(ni • v-i)n\ 

r-y 

^ 2mfm2 / f^/^ . ^; w + J_ _ I^l. + _ 16mim2(n2 ■ ^i)(i/i ■ ^2) ^» 

+(1 ^ 2), (5.5) 
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[— 4^'^''*4^^fej] ggp = superpotcntial — in — series part of {—4h'''''^ih'^k,l} 

= ( -3('^i • ^2)ni - 3r-(«i • ^'sK + — 2 K • ^2)^2 

+(1 ^ 2). (5.6) 
Below, we first consider the superpotential part of sA]^. 

A. Superpotential part 

The integrands ^k]^ - (S - dependent parts of {-ih^^^^h'' kj + 2i}S^^-^Uh'' k.i + 4/i*'',r4/i^^,/c + 24/i^^,fe6/i^'''''l}) can 
be simphfied into a form which is independent of S. For the so-obtained S-independent group, we have to find 
particular solutions of Poisson equations whose source terms are,^ 

Jill r\r-[ ■ j r\r{ r\r[ r\r{ r\rlr'[r{ 



„5' 4' 2' „7 ''I'l' 6 

r-^ r-^ 



1 1 ra 1 1 rf 1 1 1 1 



r-[r2 r-|^r2 r-j^^ra ri r[ r^rj^ r^r2 rJ-Tj rjr2 rjr2 rir2 

? 7 7 7 V 7 2 V 7 7 7 ? 7 7 7 2 7 7 7 7 7 7 7 7 7 7 

rj^r^ rj^r^ r^r^rj r^r^^ r^r^^ rjr^^ r\r\r2 rir\ r\r\ r{r\r2 r\ri r\r-Y 

rjrg rjr2 rjrg rjr2 ^'£'''2 rjrg rjr2 rir2 

r5^r2 r\r2 r\r2r2 r\r2 r\r2 r\r2r2 r\r2 r\r2 r\r{rir 



rfr2 rf rfr2 rfr2 rf rf r\r2 rjrj 



(5.7) 



It should be understood that there are Poisson equations with the same sources but with (1 ^ 2) to be solved. 

Our method to derive the superpotentials is heuristic; there are few guidelines available to find the required super- 
potentials. We proceed as follows. First, we convert all the tensorial sources into scalars with spatial derivatives. For 

example. 



r\r{rir2 

™5 ^3 



\d,^djd,Hd,i (-\ + ]-{5'^d,u + 5^dj + 5^''d,i)d,i f — ) . (5.8) 
3 ^1 2i ^1 ^2 \^^2/ 3 ^1 "^1^ \r1r2J 



(Here and henceforth, it should be understood that in general, "scalars" can have tensorial indexes carried by va and 
ri2, but do not have those by r^.) 

Second, we find the particular solutions for Poisson equations with the scalars as sources using a formula 
A{f{x)g{x)) = g{x)Af{x) + 2V/(if) • Vg(a;) + f{x)Ag{x) valid in N/B. We also use superpotential chains such 
as 

f{-3-2) ^ g^(-5-2) 
i A22 i A22 

2/(-3-4) ^ 12/(-5-4), 

where 

= ^ In 

and /(™'") satisfies A/^™'") = r5"r2 . For example, for Eq. (5.8), it is easy to find a particular solution. 



rlr{r^ri2 ^ ^ 
rl4 



^ dz\dz{dz'ldz2 3 



^What particular solutions are required depends on how we simplify the integrands. 
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where f^^' is given in [33]. Another example is 

1 



-g— = AAii^. 
rjr2 12 



■3,1) 



Following the method described above, we could find all the required particular solutions other than 

J r\r{ r\r{ rlr^ I 



In Appendix G, we list some of the particular solutions that we extensively used to derive the 3PN gravitational 
field. Other useful superpotentials arc given in [28,33,40]. The necessary particular solutions can be obtained by 
taking derivatives of these superpotentials with respect to or or Z2 or some combinations of them. 

For example, a Poisson integral of the superpotential part of — 4/i'^''*4/i'^fe,i (Eq. (5.5)) can be evaluated as follows. 
We find in N/B 



= A 



4m? 



SP 

1 



15rf 



3 ''I 



Inri 



n 



- 2m\vlv'l ( Sifclnn + 



+m\m2rl2V2 



92H-2,-3) 2 



dz\dz^ 

16mim2(ui • V2)vx-—-jr + (1^2) 



''12 ^ ^ 



dz\z2 



= ASP{t,x). 
Then the surface integral in Eq. (3.31) gives 



(5.9) 



1 



-47r 



dSk 



d(N/B) 



' ^ 5P(r,y)-5P(r,y)A 



'^'^l.M <ik> 
5rl '''''' 

l^>'"j''i''i _ 
'■ieRiri 
+ (1^2). 



|x — 7/1 dy'' 

^ — 3 In ei?i 
5 

3eiiiriri2 



1 



rl2r2 



x-y\ 

V2 + 2(ni2 • V2)n\2 + 8(fii2 • V2)n\2 In 



( 111. 
[eR2 



(5.10) 



We shall ignore the last two terms in the second to last line of the above equation, 16mf ?;}/(3ei?iri) — 
8rnlm2V2/ {'SeRiriri2) (and the same terms with the label of the star exchanged, hidden in (1 ^ 2)), because they 
have negative powers of Ra- We combine the above results and evaluate the Poisson integral of [— 4/i'^'''4/i^fe,i] gp as 



/ 

J N/B 



\x-y\ 



-ih^^''^ih'^ k,i 



SP 



12mf 
~5rf 



k„<ik> 



-V, n 



In 



{eRi) 



Ami 



vf ((ni • vi)n\ - v\) 



2 2 

+m{m2 ( — o— (m • V2)n\ H (m2 • V2)n\ + 

^ rir2 riri2r2 



riri2r2 



(ni • V2)n\2 



+ 



rl2r2 



(ni2 • V2)n 



12 



41n( ^ 
r2 



4 In 



{eR^ 



rlr2 



r2 
^1^12 



16mim2 ,^ ^ , 
+ a (^1 -^2) 



\ri2 S 



((ni2 • vi) + (n2 • vi)) 



+ n 



-^("12 • vi) - ~(^i2 • vi) - -^(n2 • -Ci) ) ) + (1 ^ 2). 



(5.11) 



Notice that in the above equation there appear ln(ei?i) and ln(ei?2). 
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In general, with the help of the superpotentials and using Eq. (3.31), we can integrate the superpotential part of 
gA]^. The explicit result is too long to write down here. ^ 



B. Superpotential-in-series part 



We could not find particular solutions for the following sources; 



' 1' 1 



rfr2 



+ A 



' 1'2 

^2 = 2^^?^- 



f (-4,-1) 

4 

ijy y;(-2 -3) 



1 



rlr2 



For these innocent looking sources, we did not find particular solutions valid throughout N/B in closed forms. Instead 
we looked for those valid near the star, say the star 1; the field we need when we evaluate the evolution equation for 
PIq and the equation of motion for the star 1 is the field around the star 1. 

Now, all the integrands classified into the superpotential-in-series part at 3PN order are found to have the following 
form (neglecting the m^, va, and f\2 dependence appearing in actual applications of the following formulas): 



d,.d.g{x)=d,.d^ 



(lnri)P(lnr2)'' 



1'2 



(5.12) 



where a and b are integers and p = 0, 1, (j = 0, 1. A, A' = 1,2. 
Then, we take spatial derivatives out of the Poisson integral, 

/ 7^9^' 9{y) = d^, I d'yji^ 



+ 9, 



dB, 



IN/ 1 

'\x-y\ 



OBa 



9{y) 

_ A' 

\x-y\ 



(5.13) 



For the remaining volume integral, we change the integration variable y into yi, namely, 2/2 = ri2 We also change 
the integration region N/B into N\/B, where A''i = — ^i| < TZ/e}, 



J N/B \x-y\ 



Ni/B 
1 



Zt Zk Zi 

21 1 11 

1 



,3 .9(2/1 + Zl) fc 

d yi— zrr - z^ 

n - yi\ 

dSkdyidym 



dSk 



ON 



dN F-2/1 

gjy) 
\x-y\ 



2! 



4! 



--^t^-^r-z? 4> dSkdyidy^dy^ 

'dN 



gjy) 
\x-y\ 



dSkdy 



ON 



gjy) 

x-y\ 



(5.14) 



where dyi = d/dy^. The surface integrals and terms expressed as • • • arise due to the difference between TV and Ni. 
See, e.g., [55]. Note that rA = x — za^ where x is the field point, while yU = y~ za, where y is the integral variable. 
For the first volume integral on the right-hand side of the above equation, we expand \ fi—yi\ as 



1 



\ri - 2/1 1 



c=0 



ri ■ yi 
riyi 



(5.15) 



''The number of terms are ~ 10^, depending on how we simplify the result. 
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where r> = max(ri,yi), r< = min(ri,yi), and Pc{x) is a Lcgcndre function of order c. 

Now we split the integration region into four parts according to where the radial variable yi is as follows: region I, 
yi G [ei?i,ri]; region II, yi G [ri, ri2-ei?2]; region III, yi G [ri2-ei?2, ?'i2+e-R2]; and region IV, yi G [ri2+ei?2,7?./e]. In 
the third integral region, angular integration is incomplete due to the body zone 2 (-B2), which the Poisson integration 
over N/ B does not cover. 

Then first for region I, 

d^y^ o V- 1 r . (In J/1 )^ 



Pc(cOS0)Pc(cOS7) , 9 9 „ 

dcos9 ^ ,99 \. m./9 11^^12 + 2/? - 2ri2yi cose* 

2«(ri2 +2/? - 2ri2yiC0s6')V2^ 



c=0 1 12 JeRl/ri2 

^ yi dtPeW (i ln(l - 2Ci + C') + Inria)' 



(1 - 2Ci + C^)^/2 

where COS7 = -(fi • ri2)Ai/ri2, t = cos9 = -{yi ■ fi2)/yi/ri2, and C = yi/ri2. 
Next for region II, 



(5.16) 



^^5(2/) = 27r > / dyi 



// \ri-yi\ - ^ "^ri ' 2/?+' ' 



Pc(cOs6l)Pc(cOS7) /w 2 , 2 o 

^ V12 + 2/1 ^ 2ri2yi cos 6^)"/^ 



^ r^Pe(cos7) rfC(lnC + lnri2)P 

"".Z^ ra+b+c-2 / , Aa+c-1 
c=0 '12 Jri/ri2 S 

1 dtPcft) ln(l - 2Ct + C') + lnri2)' 
.1 (l-2a + C')''/2 



(5.17) 



Third for region IV, 



IV Fi - Vi 



c=0 ''ri2+eR2 2/1 



Pc(c0S^^)Pc(c0S7) , o 9 „ 

^^" ^n^?2+ 2/?- 2rr2,,cos.)V2 0-(-^^ + " ^^-^^ -^))^ 



^ r^Pe(cos7) /■^/(^'-i^) rfC(lnC + lnri2)P 

'^Z^ a+b+c-2 ra+c-1 
c=0 12 •'l+«rt2/ri2 '5 

^ yi dtPe(i) (i ln(l - 2Ct + C^) + lnri2)' 



(5.18) 



(1 - 2Ct + Cfi^ 

Now for region III, the angular deficit due to the body zone 2 is determined by 

(ei?2)' =yl+ - 2ri22/i cos^o- (5-19) 

It is convenient to redefine ^ as ^ = yi/ri2 — 1. Then 6 ranges from —1 < cos^ < cos^o = 1 — C({C,), where 
a(C) = (ei?2/ri2 - C)(ei?2/ri2 + C)/2/(l + C)- Thence 
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Jiii\ri-yi\ ^ 



ri2—eR2 Hi 



r""^" , n Pc(cOs6l)Pc(cOS7) , 2 2 „ 

7_i 2'2(rf2 + 2/r - 27-122/1 cos 6')''/2 

_ ^r|P^(cos7) r 

-'^^2^ „a+6+c-2 / 

""/i "'''(2V2C-2(1 + C)/+C2)V2' ^^-2°^ 

Then summing up the above results, we obtain the following formula, by which we evaluate the first volume integral 
on the right-hand side of Eq. (5.14): 



r^Pc(cos7) rfC(ln(l + C) + lnri2)P 

— a+b+c-2 j Q , AXa+c-1 

— 12 J-eR2/ri2 y-^^) 

dtPS) ln(2 + 2C - 2(1 + Qt + C^) + lnri2)' 



,3„. 9{yi + zi) 



Jm/B \ri-yi\ 

fri/ri2 



c=0 1 '12 J<iRi/ri2 

dtPejt) (I ln(l - 2Ct + C^) + Inria)^ 

""7-1 (i-2a+c^)V2 

^ r^P,(cos7) /■^-^^^/''^^ dC(lnC + lnri2r /-^ diPe(i) ln(l - 2^ + C') + lnri2)' 

2^ „a+&+c-2 / /-a+c-l / ' 

c=0 12 Jri/ri2 S ■' — 1 

,_^ r^P.(cos7) r 

+ ^^Z^ ,a+6+c-2 / 

n '19 ^ — < 



ir./n. 7-1 (l-2a + C2)V2 

r^Pc(cos7) /■'^«2/ri2 dC(ln(l + C) + lnri2)P 

' -eR2/ri2 



(i + C)"+^-i 



d<Pc(^) (^ ln(2 + 2C - 2(1 + Qt + C') + lnri2)' 
.1 (2 + 2C-2(l + C)t + e)''/2 

'fPe(cos7) /■^/(^''i^) dC(lnC + lnri2)P 

1^ a+6+c-2 / (-o+c-l 



, „ v-l1Pc(cos7) /■' 



(1 - 2a + C^)^/2 

As an example, when the source term is r\r2/rf/r2, we have 



where 



_^ P4(C0S7) 



-^[1 2]' satisfies 



(5.21) 



/ = + 0(ei?A), (5.22) 



''''1 ijv/B -47r|f - ^ yfy2 Jn,/b -47r|fi - yi| yfy2 {^11 J 

ri2 ri2 V''i2/ Srfo V 3 V^ia 
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Thus, in the neighborhood of Zi, F^^ 2]' is the required solution of the Poisson equation in the sense of Eq. (5.24). 



In general, F^^^^ denotes a function which satisfies 



AF, 



(m,n) 
[A,c] 



'1 '2 



o 



as rA 



0. 



An appropriate value of the index c depends on how many times we should take derivatives of F^^^^ to derive an 
equation of motion. 

Now, to illustrate our method, let us evaluate the Poisson integral of [— 4/i'^'''4/i^fc,i] ggp defined by Eq. (5.6), 



52 



In ri 



92 



dzlz'[ \rlr2 



Then we evaluate the Poisson integral around z\ as 



-f 



(-2-3) 



(-4,-1) 



+ (1 



dz'[z2 
2). 



N/B -47r|f-^ 



\_—ih^^'^ih'^k,i 



= m\m2V2 



Ani\ni2 



'12 

4m 17712 



' 12 



94 [1-21 



SSP 

92 



■F, 



dz\z\ [1'21 
92 



(-2,-1) 



o f(-2,-l) 



F, 



dz\Z2 '"^'2 



(-1.-2: 



'f2 



T-2 
fi?2 



Ve^i 



ri 



1 

''2 
1 

ri 



4r2 
4r2 



-Rix) 



rfr2 



(5.25) 



(5.26) 



where R{x) is the remainder. 2] "^^ satisfies 



(ln,-3) 



In ri 



[1,2] 



Similar equations hold for F^^ 2']''"'' ^.nd -^[i 2^ 

We note that for the superpotential-in-series part, there is no need to add terms corresponding to the surface 
integral terms in Eq. (3.31). We note also that the surface integrals in Eqs. (5.13) and (5.14) do contribute to the 
field in the neighborhood of the star. 

We could evaluate the Poisson integral of the superpotential-in-series part of in the neighborhood of the star 
1 by means of the method described in this subsection. 



C. Direct-integration part 



For the direct-integration part (e.g., Eq. (5.4)), we evaluate the surface integral in Eq. (4.7) directly, while we give 
up deriving the corresponding contributions to the field valid throughout N/B in a closed form. In this subsection, 
we consider only the effect of the direct- integration part of g/i^* on the evolution equation for 

Let us define the "DIP" field g/iBip. 

g/iDip = (direct - integration part of g/i"^') = 4 / ,^ ^ sAo', (5-27) 

J N/B \x-y\ 

with 167r8A5* = (S - dependent parts of {-4/i'''''4/i^/c,; + 24/i*(''''^4/i^fc,j +4/l'^T4/l^^,fe + 24/i^^,fe6/i^'''''^})- Then in the 
derivation of the evolution equation for the direct-integration part appears as (see Eq. (F3)) 
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dR 



le 



dSk 
dSk m2r2 

Then it is sufficient to compute the following integral: 

dSk TO2''2 



T[l,k] 



27r 



■[i,fe] 

DIP 



dSk — 

IT 



2 m2r2 



"dip "t 



J3„ 8 



N/B \X 



y\ 



l]r 



d(N/B) 



V\ 



(5.28) 



(5.29) 



Straightforward calculation shows that sAg^{y) 
integral over dN. On the other hand, 



as y — > oo, and thus no contribution arises from the surface 



dSihs^ 



[fcsAg _ 4:17111712 ( mi [fe i] 2mi ;] 



dB 



V\ 



' 12 



ri 



3ri 



+ (1^2). 



(5.30) 



Therefore, the second surface integral in the square brackets in Eq. (5.29) gives no contribution to the evolution 
equation for P^q. The first integral requires special treatment, which we shall explain below. 
Now let us consider an integral which has a form 



dSk 



2 ^Ra 



disc / d yi-rz; 



fin) 



N/B 



\ri - yi\ 



(5.31) 



Here f{x) carries tensorial indexes in general, but we do not write them explicitly for notational simplicity. We call 
this type of integral the companion star integral. The first integral in Eq. (5.29) is a companion star integral with 

f{yi) = (2m2/7r)8Ag''''^^(yi + zi) (^2 must be replaced by ri2+yi in f{yi))- In the above equation, we defined disc,;/}^, 
which means to discard all the ei?A-dependent terms other than logarithms of cRa- Thus, for example. 



disc 



eRi ( ei?2 ) ''I ( £-^2 ) eRi ri 



n n \eR2 



The symbol discefl^ is introduced in Eq. (5.31) to clarify that we discard cRa dependence in the field before we 
evaluate the surface integrals in the general form of the 3PN equation of motion. 
To evaluate a companion star integral, we first exchange the order of integration. 



dSk^ disc 

rS eRA 



d^'y 



N/B 



fivi) 
\ri - yi\ 



lim disc / d 2/i/(yi) <J> dSkjzr- 

J^eiJi cRa Jj^ib JdB{ 1 2/1 



r-12 



(5.32) 



where we defined a sphere B[ whose center is z\ and radius is which is a constant slightly larger than eR\ for any 

(small) e (ei?i < r[ « ri2). 

We mention here that the procedure of exchanging the order of integrations here is motivated by the works of 
Blanchet and Faye [39-41]. 

The reason we introduced r[ is as follows. Suppose that we treat an integrand for which the superpotential 
is available. By calculating the Poisson integral, we have a piece of field corresponding to the integrand. The 
piece generally depends on eRA, however we reasonably discard such i?yi-dependent terms (other than logarithmic 
dependence) as explained in Sec. HIE. Using so-obtained iJ^-independent field, we evaluate the surface integrals 
in the general form of the 3PN equation of motion by discarding the eRA dependence emerging from the surface 
integrals, and obtain an equation of motion. Thus discarding-ei?^^ procedure must be employed at each time, when 
the field is derived and then when an equation of motion is derived, not in one time. Thus r[ was introduced to 
distinguish two species of eRA dependence and to discard cRa dependence in the right order. We show here a simple 
example. Let us consider the following integral: 



dSk- 



dBi 



N/B 



d^y 1 
^ - y I 2/1 ' 



(5.33) 
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Using A Inn = ^/rf, we can integrate the Poisson integral and obtain the "field", 

JN/B\x-y\yt V^A/ n \ J 

Since the "body zone contribution" must have an eRi dependence hidden in the "moments" as 47rei?i/ri(+0((ei?yi)^)) 
(see Sec. Ill E), the terms — 47rei?i/ri+0((ei?A)^) should be discarded before we evaluate the "equation of motion" (the 
surface integral in Eq. (5.33)). The surface integral gives the "equation of motion", 



eRi 

On the other hand, we can derive the "equation of motion" by first evaluating the surface integral over dB[ 



/ dsAl ^/ ^^^^ 

/asi n Jn/b \x-y\vi Jn/b Vt JaB[ 



fi \ri - 2/1 1 



= 167r^ 









Jr{ yi 



lQ^^l^(^^l)+l-^^). (5.35) 



Thus, if we take dBi as the integral region instead of dB[ in the first equality in Eq. (5.35), or if we take limrj_>eHi 
without employing disCeij^ beforehand, we will obtain an incorrect result. 



IStt^ In 

which disagrees with Eq. (5.34) 



(eii)' 



Now let us return to Eq. (5.32). To evaluate the surface integral, we expand the integrand, supposing r'l is small, ^ 

1 



79. 



fdB[ ^^'^ Kni - yi I \ri2 + r[ni 



'a+2 ( 'b ,,fc 
'1 Vl 



a=0 
6=0 



. - I ■ ,,,'6+1 

12 [Vl ' I 



-And,, Yi-ir ~ ^}'\ \, n<i^^>N<^-> J ^^,r;V \ , (5-36) 



where Nl = yl/yi, and in {f,g} in the above equation, / denotes the result for r[ < yi and g denotes the result for 
r[ > y\. In the last equality of Eq. (5.36), we used the following formula [58]: 



1 

47r 



-N^ 5a,b+l+ , TVi ^l "12 <5a+l,f,. (5.37) 



(2a +1)!! ^ "'"^^ ' (26+1)! 

Substituting Eq. (5.36) into Eq. (5.32), we establish a formula 

^0!! = (-1)!! = 1. 
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dSk disc / d y 



fim) 

dBi »'2 "■"■■4 Jn/b " "In-yil 

= lim disc47ry(-ir i^"~y, 

H^eflA eflA ^ (2a + 3)a! 

0=0 ^ ' 

'2a+3 



<kMa> 



\ '12 

'12 

^12 



(5.38) 



where = UB2. 

Now since wc expect from the discussion in Sec. HIE that an equation of motion does not depend on sRa and hence 
r'l, in the formally infinite series in Eq. (5.38), only finite terms for which the volume integral makes r'^ dependence 
disappear do contribute. The r[ dependence of the above integral can be examined by the behavior of the integrand 
near dB[. Thus, to examine r[ dependence, we expand f{yi) in the small yi,^ 



P=Po ^1 p 



(5.39) 



where po, the lowest bound in the sum, is a finite integer. Then, r[ dependence of the four volume integrals in Eq. 
(5.38) can be evaluated as 



'a-p+4 



disc / dy,fim)^ = E /(^i • ^12) ' -^ Si-P^a - r[^'^'h-p,a Inr; 



eRA 



„'2a+3 



disc / dyifiyi)-^-^ 



E /(^i • "12) I ^^p_/ 2-p^a - r^+'<52-p,a Inr'A , 



P=Po 



'a—p+i 



disc / dyifiM^' = E /(^l • "12) ' ^/ p-4#a + .5p_4,a In ^ 

' a—p+5 



P=Po 



disc / dy^f{y^)T[^y1+^ = ^ f{N^ ■ n,^) ' ^J p-3^a + r^H^.^^a In 
•^^^ JeRi rzt^ P \a-p + 3 v^-Ki 



P=PO 



(5.40a) 
(5.40b) 
(5.40c) 
(5.40d) 



where fia^h = for a = 6 and 1 otherwise. Note that a > 0. Since we take cRa and discard cRa dependence, in 
the above four equations the only terms which we should retain are 



f ^'2a+3 ^1 

disc lim disc / dyif{yi) ^ = } -Sa,p-i f{Ni ■ ni2), 

cRa r[^eRi cRa Jr> Wi ^ 2p - 5 p 

"'i p=po 

f /2a+3 ^ ^ ^ 

disc lim disc / dyif{yi) \ = V -Sa,p-5fiNi -7112), 

€Ra rj— »eiJi cRa J j.i y^ Zp — I p 



P=Po 



disc lim disc / dyif{yi)y^+^ = 0, 

eRA r[^eRi cRa J 



disc lim disc / dyif{yi)r^y1+^ = - V -5a,p-5 f{Ni ■ ni2). 

cRa r{^eRi eRA J^ili ±1 ^ P 



(5.41a) 

(5.41b) 

(5.41c) 
(5.41d) 



Substituting Eqs. (5.41a), (5.41b), (5.41c), and (5.41d) into Eq. (5.38), we obtain 



'gAg* has no logarithmic dependence. 
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■4 



disc 



,3 fjyi) 

IT "-iiov. I a y— 



T 

dSk disc 



isc / 



J/1 1 



a=po— 4>0 

+ 477 5^ (-1) 

a— po— 5>0 



g (2a +1)!! 
2(2a + 3)2o! 



dON, / (iVi • ni2)7V<^»>a, 

0+5 



'12 



' 12 



<fcMa> 

'-12 



(5.42) 



where we applied disce/j^ to the left-hand side of Eq. (5.38) to make it clear that we discard the ei?A-dependent terms 
(other than the In cRa dependence) when deriving an equation of motion. 

To illustrate our method, let us take [— 4/i'^''*4/i^fe,(] j-,jp defined by Eq. (5.4) as an integrand, for instance, in Eq. 
(5.42), 



2m2 



'[3 



k,l 



DIP 



where the vertical strokes denote that indexes between the strokes are excluded from (anti)symmetrization. We 
expand this integrand around z\ as 



2m2 



TT 



-4,11' Un k,i 



DIP 



kl\ 



32mlml 



^-nia^i' - 2(ni • vi)n}ln{^2 ~ 2(«i • «i2)«i^^i') 




(5.43) 



DIP 



for which po = 4, does not contribute to the 3PN evolution 



Thus we see that the integrand dy [—4/1' 

equation of P^q since the angular integration of the integrand in Eq. (5.42) gives zero. 

Returning to the evaluation of s^Idip' apply Eq. (5.42) to f{yi) — (2m2/7r)8A^'''*^' (yi + zi). Expanding 
8Ag'''^'(^i + zi) around zl, we find that po = 4. Evaluating the angular integral in Eq. (5.42) for 8Ag''''^'(y'i + zi), 
however, results in zero. 



VI. HARMONIC CONDITION 

In this section, we check a part of the harmonic condition, 

<8h^^,r + <8/l",i = 0. (6.1) 

As explained in the previous section, we split the nonretardcd part of the N/B contribution 8^]v/Bn=o ^^^'^ three 

groups: g/isp (t^® superpotential part), s/^ssp (^^^ superpotential-in-series part), and s/idip i^^^ direct-integration 
part). We could obtain an explicit form of only s^sp' body zone contribution sh^, and the second time derivative 
term in the retarded field 8'ijv/Bn=2 (*^^ ^ast term in Eq. (5.1)). Our trick is to transform the left-hand side of Eq. 
(6.1) into the following form: 

Ti{T , x) = <8^^^,T + <7h'^\i + 8h^,i + s/^SP,* + 8^JV/Bn=2,j 

J N/B F - y\ oy Jd(N/B) F ~ y\ 

where sAggp is the integrand corresponding to the superpotential-in-series part. Now not surprisingly, we could 
evaluate the integrals explicitly using superpotentials. For example, 

Ai6.8A§|p = A ( 24mf m2yl(yl • .-2) . 
dt dy^ \ yfy2 

= A (24m?m2.2^ (-^a../^-'-^) + f a../^-^-^) - ^^.Z^"^'-^))) 

+ •••• (6.3) 

In this manner, though wc could not evaluate Poisson integrals explicitly and consequently the field g/i"^* valid 
throughout N/B is not available in a closed form, we checked the harmonic condition in the sense that T-L{t,x) = 
throughout N/B. 
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VII. THE 3PN MASS-ENERGY RELATION 



As explained in Sec. IV, the direct- integration part does not play any role in the evaluation of the evolution equation 

of 

^Ae- Thus wc can take the same method as in the evaluation of the 2.5PN equation of motion. We first express 
io6]^ explicitly by substituting the field derived previously into Eqs. (F2) and (F3). Then evaluating the surface 
integrals in Eq. (4.7), we obtain the evolution equation of P^q, 



le 



dT 



<3PN 



+ € 



217111712 



' 12 



4mim2 



' 12 



[4(ni2 • vi) - 3(ni2 • V2) 



-^(ni2 • ^2)^ + ^Vi(ni2 • V2) + 6(ni2 • Ui)(ni2 • V2f 



2vl{ni2 ■ vi) + 4{vi ■ V2){ni2 ■ V) + hvl{ni2 ■ V2) - 4t;|(ni2 • vi) 



+ — (-4(ni2 • V2) + 6(ni2 • vi)) + — (-10(ni2 • vi) + ll(ni2 • V2)) 
ri2 ri2 



gmiTO2 



-vl + 2v{vi + Av^ (ni2 ■vi) + 



:ViV2 + IV2 (ni2 • V2) 



+ {2vl + Avfj (ni2 • vi){vi ■ V2) - {2vl + 8v|) (ni2 • V2){vi ■ V2) 
+ {Zvl + I2vl) (ni2 • i?i)(ni2 • Vif - {^vl + Uvl^ (ni2 • ^2)^ 
+ 2(ni2 • V2){vi ■ V2f - 6(ni2 • '?i)(ni2 • V2f'{vi ■ V2) + 6(ni2 • V2f{vi ■ V2) 



y (?^i2 • vi){ni2 ■ V2Y + y (^12 • ■^^2)^ 



mi 



ri2 



117 



— ( ( -A2vi - ^v'i 1 (ni2 • vi) + 60(ni2 • vif 



/137 2 37 2\ ^ s 

( — ^'i + Y^^2 1 (^2 -^^2) 

297 219 

(^12 • vi){vi ■ V2) ^(^12 • V2){vi ■ V2) - 151(ni2 • ■yi)^(ni2 • V2) 



4 ' ' ' ' 4 

+ 109(ni2 • wi)(ni2 • V2f - 23(ni2 • V2f) 

+ 26(ni2 • v\){v\ ■ V2) - 28(ni2 • V2){vi ■ V2) + 2(ni2 • '?i)^(ni2 • V2) 
+ 16(ni2 • vi){ni2 ■ V2f - 20(ni2 • ^2)^) 



mf / 33,^ ^ , 13,^ 

— K2 • ^^i) - y("i2 • ^2, 



miTO2 /^35 ^ 17 ^ 
^("12 ■vi) + -^{n\2 ■ V2) 



' 12 



TTln 



' 12 



23_ 



+ — -12(ni2 • v-i) + — (ni2 • V2) 



(7.1) 



where V = v\ — V2- 

Remarkably, we can integrate Eq. (7.1) functionally. 



PIq = mi (1 + e22ri + eVi + e^Fi) + O(e^) 



(7.2) 



with 



1 2 , 3m2 
2ri = 7,^1 + — , 

2 ri2 

„ 3rn2 ^ 2 , 2rn2 2 , 7rn- 
4ri = -7^— (ni2 • t;2) H ■?^2 + 



2ri2 



n2 



2ri2 



2 2 



4ms 
r-12 



^ s 3 4 7m^ 5mim2 
■(vi-V2) + -t^i + 7r:2 



2r2 



12 



(7.3) 
(7.4) 
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m\m2 2lTOim2 5to2 



o„3 
^'12 



4r32 



2rf2 16 



+- 



ml ( 45 
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1 _ 



v{ + ^t;2 + ^("12 • v\Y - 19(ui • V2) - (ni2 • ui)(ni2 • ^2) - 3(ni2 • v^f 



mim2 / 43 2 53 2 69,^ ^,3 53,^ 85,^ \ 

l^y^i + -§-^2 - -§-('^12 ■ ^1) - • ^2) + ^("12 • vi){ni2 ■ V2) 

69, — Tn'2 /^33^ *-»22 r\4 n 1 f —* —* \ A 1 / — *\ 

- -g-("'i2 • ■^^2) ) + ~ V "s"^^ ^^^^ ^2 ~ (Ui • V2) - ^V2{Vi ■ V2) 

- ^^^l(^12 • V2Y - ^Vl{ni2 ■ V2f + 2(Ul • V2f' 

9 

+ 2(ni2 • V2f{vi ■ V2) + -(ni2 • V2) 



(7.5) 



The mass-energy relation for the x part up to 3PN order is given in Appendix E. Equations (7.2) and (E2) give the 
3PN order mass-energy relation. 

1. Meaning of PXe 

In this section, we suggest an interesting interpretation of the mass-energy relation. First of all, we expand in an e 
series a four-velocity of the star A normalized as g^vu'^u'^ = — e~^, where u\ = u\v\. We have 



4" 4" " 16 16 128^ 'a 



An k 



1 115 



6 



16 



(7.6) 



This is a formal scries since the field should be evaluated somehow at za while the metric derived via the point-particle 

description diverges at za- 

Now let us regularize this equation with Hadamard's partie finie rcgularization (see, e.g., [41] in the literature of 
the post-Newtonian approximation). Evaluating with this procedure Eq. (7.6) and y/—g expanded in e up to O(e^), 
then comparing the result with Eq. (7.2), we find at least up to 3PN order 



Ple = rnA[V^u:,]T. 



(7.7) 



In the above equation, [/]^* means that we regularize the quantity / at the star A by Hadamard's partie finie or 
whatever regularizations which give the same result. 

We emphasize that we have never assumed this "natural" relation in advance. This relation Eq. (7.2) has been 
derived by solving the evolution equation for PJq functionally. 



VIII. THE 3PN MOMENTUM- VELOCITY RELATION 



We now derive the 3PN momentum-velocity relation by calculating the Q\ integral at 3PN order. From the 
definition of the Q\ integral, Eq. (3.24), 



Q\ = / dSk (loA^v'^ - t^lioA]7) + O(e^), 

J dBA 
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we find that the calculation required is almost the same as that in the equation for dP^Q/dr. Namely, s^dip does 
not contribute to the Q\ integral due to (i) the antisymmetry of the direct-integration part (see Eq. (5.28)), (ii) the 
behavior of the direct-integration part at infinity (see Eq. (5.29)), and (iii) its behavior around za (see Eq. (5.42). 
Note that s^y^iT, za +2m)2/a ~ ^ly\ neighborhood of the star A). Therefore, we need to compute the surface 

integrals in the definition of Q\ using the field up to 2PN order, ^h7j^, s^^N/Bn=2^ s^isP' ^'^'^ s'lssp- show only 
Q'aq liere, 

<6Qie - -e e - y-^r,, J _ -e - (^-m,a, J , (8.1) 

where it should be understood that a\ in the last expression is evaluated with the Newtonian acceleration. We show 
Q\-^ in Appendix E. 

In the field, Q\ of O(e^) appears at 4PN or higher order. Thus up to 3PN order, qQ\ affects the equation of motion 
only through the 3PN momentum-velocity relation. As stated in Sec. IV, the nontrivial momentum-velocity relation 
which affects the equation of motion is that of the Q part. This in turn motivates us to define the representative 
point of the star A, z\, by choosing the value of D\q. We do not take into account £)\^ in the definition of z\. 

Now with Q\q in hand, we obtain the momentum-velocity relation, 

Pl,,P;^.-^±(i^l^W^. (8.2) 

This equation suggests that we choose 

D\&{r) = ^'\<a\ = e'S^Aeir). (8.3) 

For a while henceforth, we shall define z\ by this equation. We shall later give a more convenient definition of z\. 

Finally, we note that the nonzero dipole moment D\q of order e** affects the 3PN field and the 3PN equation of 
motion in essentially the same manner as the Newtonian dipole moment affects the Newtonian field and the Newtonian 
equation of motion. From Eqs. (3.13), (3.14), and (3.15) (or Eqs. (AS), (A9), and (AlO) for more explicit expressions), 
we see that 5\q appears only at wK^^ as 

/^^n..e=4ei° E %^ + 0(e"). (8.4) 

Then the corresponding acceleration becomes 

"^laiU^e - -e — 3 "12 +e — 3 «i2 -^—^ZT- (8.5j 

O-T 

The last term comes from the momentum-velocity relation Eq. (8.2) and compensates the Q\q integral contribution 
also appearing through Eq. (8.2). 
Note that this change of the acceleration does not affect the conservation of the binary orbital energy, 



^AS^-VA-j^ 



(8.6) 



IX. THE 3PN GRAVITATIONAL FIELD: N/B INTEGRALS 

To derive a 3PN equation of motion, we have to have lo/i^^ + &h^k besides s/i^' and the 2.5PN field. The 2.5PN 
field is well known [28]. (See paper II for ^h^^ .) At 3PN order, while the body zone contribution <iah^g + <shgk is 
easily found as Eqs. (AS) and (AlO) with the help of the 3PN mass-energy relation, Eqs. (7.2) and (E2), it is quite 
difficult to evaluate the Poisson-type integrals over the N/B region. The problem is again to find superpotentials 
required in the derivation of the field. Thus we have not evaluated all the Poisson integrals and we have applied the 
method used in the evaluation of the s^^' contribution to dP^/dr. 

We shall deal with the nonretarded field from the next subsection, 

4 / jf^ {,oA^J{T,y)+sA%k{r,y)) • (9.1) 

J N/B F ~ y\ 

In the last subsection, we consider the second and the fourth time derivative terms in the retarded field. 
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A. Superpotential part 



Using the 2PN field, we first write down iqA]^ + s^%k explicitly, and simplify it to remove its S-dependence as 
much as possible. Then as in Sec. V, wc split the integrand into two groups: the S-indepcndent group and the 
S-dependent group. For the S-independent group, we could find the superpotentials required other than those which 
have the following sources in the Poisson equations: 



1 1 1 1 r2 1 r{r{ r\r{ r2r\r{ r\r{ r\r{ r2r{r{ r\r{ 



I r°r2 rjrf r\r2 r° rfrj rfr^ rf rfr2 r°r2 rjr^ rfrj rfrj rfr^ 

r2r\r{ rir\r{ rfr^^r^ r\r\r-[ r\r2 r\r2 rfr^^rj r\r2 rlr^rj r2r\r2 1 

' ' 3 ' 5 ' ~6 ' 5 2 ' 5 ' 4 3 ' 3 ' ( • 

ri r2 r"2 r°r2 r°r2 rjrl rf ri I 



(9.2) 



It should be understood that there are Poisson equations with the same sources but with (1^2). We shall treat the 
integrand corresponding to the list (9.2), the superpotential-in-series part of iqA]^ + sA^fc, in the next subsection. 

For the remaining integrand of loAy + sA^fc (the superpotential part of iqA]^ + gA^fc), we need to find the 
superpotentials whose sources are 

11 1 1 r2 r2 r| 1 1 ''i 



^1 '^i ^1^2 ^1^2 Ti 7^ rir2 r-^r2 
rlrl rlrl r\r{ r2r\'r{ r2r\r{ r\r{ r\r{ r\r{ r2r\r{ r2r|r^ r\r{ r\r{ r\r{ r\r[ 

7 4 ) 7 2 ' ~7 ' 7 ' 7 ' ~fi 5 ' ~5 4 ' ~5 ' 5 ' 5 ' ~A 5 ' ~3 5 ' 5 ' 3 ' 

r{r| rjrj r{r2 r{ r{ r°r2 rj^r| rj^r2 rj rjrj rjr2 riTj rir^ 

r}r^ rir\r{ rir^rj r\r\r{ r\r\r{ r\r!2 '*i''2 fl^l^i ^2''^i^2 '''1^2 ^i''^2 ^1^2 ^2r\r2 



K ' 4~' 4 ' 4 ' ~3 ' 3 ' 2~' 2 • 
j,o ^4 ^4 ^Jj,^ 

^1 ^1 ^^^1 ^1 ^1 ^1 ^1 '^1 ^1 ^2 '^2 '^1 '^1 ^2 '^2 1 



1'2' 



(9.3) 



Here we did not list the sources for which we already find the superpotentials (see the list (5.7)). 

We could derive the superpotentials corresponding to the list (9.3) using the procedure described in Sec. VA. 
Useful particular solutions are given in [28,33,40] and in Appendix G. Those superpotentials enable us to calculate 
the Poisson integral with the superpotential part as the integrand. We cannot write the result down here because of 
its enormous length. 



B. Superpotential-in-series part 

The superpotentials having the sources listed in (9.2) could not be found. Thus we employed the method explained 
in Sec. VB. First we transform tensorial sources into scalars. For example. 



rlriri 1 9^ 



3 dzizi r2 



A 



1 

63 dz\zi 



214 



rl2f^''-'^ 



3rl 
7r2 



3r?2/ 



K5.-3) 



K7-5) 



We apply Eqs. (5.13), (5.14), and (5.21) to the sources in (9.2). The (scalar)sources derived by making the sources 
in the list (9.2) scalars to which we apply Eq. (5.21) are; 
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(9.4) 



and the same functions with their labels 1 and 2 exchanged. For these soiirces, we could evaluate the Poisson integrals 
in a similar sense to Eq. (5.24), and as a result we obtain in the neighborhood of the star 1 the field corresponding 
to the superpotential-in-series part. 
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C. Direct-integration part 



We now consider the DIP field contribution to an equation of motion. At 3PN order, the DIP field appears in the 
integrands of the surface integrals of the general form of the 3PN equation of motion as (see Eqs. (F3) and (F4)) 
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(9.5) 
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Here we added the s/idip contribution. (Note that in Sec. VC, we evaluated the s^dip contribution to the evolution 
equation for PJ©, but not to an equation of motion.) Then the DIP field contribution to a 3PN acceleration denoted 
by a^Dip becomes 
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where we used a relation 
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(9.7) 



All the integrals in Eq. (9.7) other than the first four terms can be easily evaluated. We now explain how to evaluate 
the first four integrals. 



1. Main star integral 

For the first integral in Eq. (9.7), wc change the integration variable y to yi and also change the integration region 
A'' to A^i using Eq. (5.14). The surface integrals over dN can be easily evaluated. Note that 2/2 in the integrands 
must be replaced by fi2 + yi ■ For the remaining volume integral, we use computationally the same method as the 



26 



one employed by Blanchet and Faye in [40,41]. Let us consider the following integral, which we call the main star 
integral: 



1 

47r 



dBi 



dfl disc 

eRa 



d^yi 



fivi) 



m/B 47r|fi-yi|' 



(9.8) 



(For the definition of disCeij^, see Sec. VC.) With the notice given below Eq. (5.32) in mind, we first exchange the 
order of integration, 
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^ for n > yi, 
^ for n < yi. 
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/asi \ri-yi\ 

Next we make a symmetric trace-free decomposition (STF decomposition) of the integrand on the indexes of ni , 

/(yi) = l^ff;(cos^,t/i)nf''>, (9.11) 
;=o 

where cos 9 = — ni2 • ni. Note that gi{cos9,yi) is not necessarily a scalar. In general, gi{cos6,y\) is a tensor whose 
indexes are carried by va, ri2, or some combinations of them. Then substituting back the STF-decomposed integrand 
into Eq. (9.9), we have 
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where (below, n and N are unit vectors) 
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J Att J 47r 



(9.12) 



(9.13) 



was used. Notice that in Eq. (9.12), when yi e [ri2 — ei?2, ru + ei?2] the body zone 2 prevents the angular integration 
region from being complete. The angular deficit is given by Eq. (5.19). 
Now let us give an example. Take d^i{\n S)dj{l/ri) as an integrand, 
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(9.14) 



(9.15) 
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The integrand (In S)dj{l/ri) was taken in [40] as an example of Hadamard's partie finie of Poisson integrals. The 
results are the same. A difference in the two results here is simply in definitions of the lower bounds of the integration. 
Here in our case we set cRa as the lower bounds. 

Now we return to the evaluation of the 3PN gravitational field. Applying Eq. (9.12) to the first integral in Eq. 
(9.7), where 



/(yi) 
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we could evaluate the main star integral contribution to a 3PN acceleration. 



2. Companion star integral 

The second, the third, and the fourth terms in Eq. (9.7) have a form of a companion star integral and thus can be 
evaluated by the method described in Sec. VC. The integrands have no logarithmic dependence and thus admit an 
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expansion in the form of Eq. (5.39), where we found po = 5. Then using Eq. (5.42), we obtain the companion star 
integral contribution, 
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D. Retarded field 

At 3PN order, we have to evaluate the integral over N/B in the second retardation expansion term, 

2^ f d^y\x-y\{sAh^ + eA%k). 

OT Jn/B 

This integral can be evaluated via the super-superpotential f{y) satisfying gA^y + eA^fc = A^/(y) as 
Jn/B 

= -8^/(x) 

+ / dS, \\x - y\d,Af{y) - A/(y) + ^—d^fiy) + ^ • 
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The superpotential of sA]^ + eA^fe is easily found. (In the derivation of s^i^^ and g/i*-' wc found the superpotentials 
required here.) The super-superpotentials we could not find are only those with l/(r^r2), l/(rir|), ri/r2, and r2/rf 
as sources. The corresponding integrands are 

f d^y f 21m\m2 15m^m22/2 21mim\ lbmim\yi\ 

Jn, 



dr'^ Jn/B ^t^\x - y\ 



yly^ 



'^rl^yl yiyl "^rl^yl J 



Thus we use the method explained in Sec. VB and obtain the field near the star 1, which is sufficient to derive the 

equation of motion. 

The N/B integral appearing in the fourth retardation expansion term can be evaluated via the following super- 
super-superpotentials: 
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It is straightforward to evaluate the surface integrals. 



\x-y\ 



\x- y|3 



f{y) 



29 



X. THE 3PN EQUATION OF MOTION WITH LOGARITHMIC TERMS 



To obtain a 3PN equation of motion, we evaluate the surface integrals in the general form of the 3PN equation 
of motion Eq. (4.8) using the field s^^^, <6^'"'j the 3PN body zone contributions, and the 3PN N /B contributions 
corresponding to the superpotential part and the superpotential-in-series part. We then combine the result with 
the contribution from the direct-integration part. For a computational check, we have used the direct-integration 
method (the method with which we evaluate the direct-integration part) to evaluate the contributions to the equation 
of motion from all of the 3PN N/B nonretarded field, 8^]v/Bn=o ^'^'^ '^o^T'n^/ Bn=o + 8^jv/Bra=o'' assuming that 
they had belonged completely to a direct-integration part. As expected, we obtained the same result from two 
methods: the direct-integration method and the direct-integration method plus the superpotential method plus the 
superpotential-in-series method. 

In the evaluation of the body zone field h'^ (shown as Eqs. (A8), (A9), and (AlO)), besides the explicitly seen 
energy monopole terms in these fields which must be converted into mass monopole terms via the 3PN mass-energy 
relation, the effects of the Q^'* and i?^'*"* integrals appearing in the 3PN field are properly taken into account through 
Eq. (A4). eQke given in Eq. (8.1) affects a 3PN equation of motion through the 3PN momentum-velocity relation. 
Since we define the representative points of the stars via Eq. (8.3), we add the corresponding acceleration given by 
Eq. (8.5). Furthermore, our choice of the representative points of the stars makes appear independently of D\q 

in the field, and hence 4D\^ affects the 3PN equation of motion. In summary, <4Q^'^, <4-R^'''' , eQ^e' ^Ae' ^^'^ '^^Ax 
contributions to the 3PN field can be written as 

lo/i"" + sh\ = 4 V 4 fee + iD'X^ + iRf - Lr'a) +■■■, (10-1) 

where • • • arc other contributions. The effect of this field on the equation of motion is properly taken into account 
via Eq. (A4) with Cqr replaced by Cqr -|- Cd^ + Ca^e: where C&^^ = 2/3 and Cd^ = -350/9. (See Eq. (8.3) 
for Cs^Q and Eq. (E3) for Cd^-) On the other hand, qQ\q and 5\q affect the equation of motion through the 
momentum- velocity relation in Eq. (4.8), 



<3PN ^'^ ^'^^ 

but cancel each other out, since we chose Eq. (8.3). We note that there is no need to take into account an effect of 
through the momentum-velocity relation of the x part on the equation of motion since we evaluate the general 
form of the 3PN equation of motion Eq. (4.8) to derive a 3PN equation of motion. Then these contributions to a 
3PN acceleration can be summarized into 

(the contribution to m\a\ from <4(5^'', <iR^^^-' ,qQ\q, 

3 2 2 3 

^'12 

_ 6 118 mfmj , .lis m\ml i 
9 9 r^^ 

Collecting these contributions mentioned from the beginning of this section, we obtain a 3PN equation of motion. 
However, we found that logarithmic terms having the arbitrary constants cRa in their arguments survive, 
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where the acceleration through 2.5PN order, {dv\/ dT)<.2.5Fn is the Damour and Deruelle 2.5PN acceleration. In our 
formalism, we have computed it in paper II. The "• • •" stands for the terms that do not include any logarithms. 
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Since this equation contains two arbitrary constants, the body zone radii Ra, at first sight its predictabihty on the 
orbital motion of the binary seems to be reduced. In the next section, we shall show that by reasonable redefinition 
of the representative points of the stars, we can remove Ra from our equation of motion. There, we show the explicit 
form of the 3PN equation of motion we obtained. 



XI. THE 3PN EQUATION OF MOTION 

The following alternative choice of the representative point of the star A removes the cRa dependence in Eq. (10.4): 

Axi /_\ ,422 3 i / ''12 \ _ Aci , ,4ri / \ _ Axi 



^Ae,New(T) = e'^S^Asir) - e^ym^a^i In (^-^ j ^ e^J^^eW + e^c^^mW ^ e*^\(r). (11.1) 

Note that this redefinition of the representative points does not affect the existence of the energy conservation, as was 
shown by Eq. (8.6). We can examine the effect of this redefinition onto the equation of motion using Eq. (8.5) (use 
5\i^ instead of 6^0)- Thence we have 
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+ 8(ni2 • Vfr^y2 - 2(ni2 • V)V') ■ (11.2) 

Comparing the above equation with Eq. (10.4), we easily conclude that the representative point z\ of the star A 
defined by 



(r) = / dS{y' - z\(t)W^(t, y^) = e^Sl^ir) (11.3) 

JBa 



obeys an equation of motion free from logarithms and hence free from any ambiguity up to 3PN order inclusively. 

We mention here that Blanchet and Faye [40] have already noticed that in their 3PN equation of motion a suit- 
able coordinate transformation removes (parts of) logarithmic dependences of arbitrary parameters corresponding 
(roughly) to our body zone radii.'' It is well known that choosing different values of dipole moments corresponds to 
the coordinate transformation. 

By adding mia\\sj^^^ to Eq. (10.4), we obtain our 3PN equation of motion for two spherical compact stars whose 
representative points are defined by Eq. (11.3), 



'^Unlike our case, their coordinate transformation does not remove the logarithmic dependences of their free parameters 
completely. The remaining logarithmic dependence was used to make their equation of motion conservative. 
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in the harmonic gauge. 

Now we Hst some features of our 3PN equation of motion. In the test-particle hmit, our 3PN equation of motion 
coincides with a geodesic equation for a test-particle in the Schwarzschild metric in the harmonic coordinate (up to 
3PN order). With the help of the formulas developed in [42], we have checked the Lorentz invariance of Eq. (11.4) 
(in the post-Newtonian pcrturbativc sense). Also, we have checked that our 3PN acceleration admits a conserved 
energy of the binary orbital motion (modulo the 2.5PN radiation reaction effect). In fact, the energy of the binary E 
associated with Eq. (11.4) is 



E = -miWi 



mim2 



2ri2 



m\m2 



+ 



mim2 
2ri2 



3vl 



'^{Vl-V2) 



-(ni2 • ui)(ni2 • V2) 



34 



5 g mf 777,2 \%m\m2 



TO?m2 / „ 2 7 2 29,^ ^,2 13/_ -N/-. -N ^ n2 

+ + -t;2 + yl'^ia • vi) - -^{rii2 ■ Vi){ni2 ■ V2) + (ni2 • 

+ ( ^(rii2 • vxf{ni2 ■ V2) + 7(ni2 • Vxf'{ni2 ■ V2f - ^(ni2 • Vi){ni2 ■ V2)vl 

13 21 13 

^(«12 • V2Yvl + yt^i + -^{ni2 ■ Vifivi ■ V2) + 3(ni2 • Ui)(ni2 • V2){Vi ■ V2) 

. 2/-> -> \ I /-> - \2 . 31 2 2 



yUl(tTl-tT2) + y(tTl-tT2) +— W1W2 



35 o 3mf?7i2 4697Tifm2 
128™^"^ + ^^fT ^ 



2 



/547_ ^.2 3115 _ ^, 123^ _ ^ , _ ^ . 
I (ni2 • vi) ;ti-("'12 • Vi){ni2 ■ V2j rr— (rti2 • fi)(ni2 • V/) 



2r?2 V 6 24 -^-.v— 32 

575 2 4l7r% ^ ^ 4429 _ 

m?r772 / 437,^ 317,^ ^ 301 

+ I ^("-12-1'l) + ^('^12 • 'yi)(ni2 • f2) + 3(7^12 • W2) + 

337,^ ^ , 5 , 

- -Y^{vi-V2) + -v', 

Tfli7n2 ( 5 , ^ \5/-t \ ^ / -t -* \4/-t -» \2 ^ / -* xS/-" -» \3 

H -T?('^12 • ^'ij («12 • U2) - T7-("12 • Vi) (ni2 • ti2) - ;^(«12 • I'l) (ni2 • V2) 

ri2 V 1" 16 32 

19 15 3 

+ Ye^"^^ ' • V2)vl + Y^(?^12 • 't^l)^(ni2 • t52)^^^? + ^("12 • Vi)(ni2 • V2fvl 

+ Y^("-12 • ^^2)^^! - Y^("^12 • Wl)("^12 • V2)vl - 2(7712 ' W2)^l'i + ^""l ~ ^^("12 ' Vi)^{vi ■ V2) 

15 

- {ni2 ■ vi)^{ni2 ■ V2){vi ■ V2) - ;^(^i2 ■ ^^1)^(^12 • V2f'{vi ■ V2) 

45 5 11 

+ Ye^^^^ ' ■ '^2) + ^(^12 • Vi){ni2 ■ V2)vl{vi ■ V2) + -j{ni2 ■ V2f'vl{vi ■ V2) 

- • ^2) - ^("-12 • ^\f{^\ ■ + ^("-12 ■ f^l)("^12 • V2){vx ■ V2f + ' ^if 

+ ^(vi • V2f - Y^(^i2 • vxfv\vl - p(«i2 • t/i)(ni2 • t72)t;^t;^ + '^■^A^l - ^'"I'^I'k^^ ■ ^2) 

m\m2 ( 49,^ _ ,4 75,^ ^ 187,^ ^ x2/^ n2 H 4 
H 2 — — S"("i2 • ^1) + "^("12 • ^1) ("12 • "^2) 3-("-i2 • I'l) ("12 • V2) + — Wi 



' 12 



2 



247 49 81 21 

+ — (ni2 • 'iTi)(ni2 • •?2)^ + -^-(^la * "i^O^^^i + -3-(^i2 • 'i^i)(ni2 • v^v\ - —(^12 • '^2)^^^? 

15 3 21 
^(^12 • • ^1) - 2(^12 • ^\){n\2 ■ V2){vi ■ V2) + -^{n\2 ■ V2f{vi ■ V2) - 27v1{vi ■ V2) 

/ -* -* \2 49 . -* \2 2 I ->■ -> \ / ^ \ 2 ^ / -* \2 2 2 2 

+ • ^'2) + ^("12 V2- ^("-12 • ^'l)(^12 • V2jV2 + ^("12 ' ^'2 ) ■*'2 + "^^1^2 



28(iTi • V2)vl + 



+ (1^2) + 0(e^). (11.5) 



This orbital energy of the binary is computed based on that found in Blanchet and Faye [40], the relation between 
their 3PN equation of motion and our result described in Sec. XII 1 below, and Eq. (8.6). (After constructing E 
given as Eq. (11.5), we have checked that our 3PN equations of motion make E conserved.) 

We note that Eq. (10.4) gives a correct geodesic equation in the test-particle limit, is Lorentz invariant, and admits 
the conserved energy. These facts can be seen by the form of ai|<5^,„, Eq. (11.2); it is zero when mi 0, is Lorentz 
invariant up to 3PN order, and is the effect of the mere redefinition of the dipole moments which does not break 
energy conservation. 
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Finally, we mention here two computational details. In the course of calculation, z\, Z2 appears independently, that 
is, not in a form as fi2. This can be seen, for instance, from Eq. (5.14). As another example, the surface integral over 
the near zone boimdary in Eq. (3.31) in general gives terms depending on za explicitly. All such za dependences, 
when collected in the equation of motion, are combined into r\2 ■ We have retained during our calculation 7^-dependent 
terms with positive powers of or logarithms of TZ. As stated below Eq. (3.9), it is a good computational check 
to show that our equation of motion does not depend on TZ physically. In fact, we found that Tl-dependent terms 
canceled each other out in the final result. There was no need to employ a gauge transformation to remove such TZ 
dependences. 



XII. COMPARISON AND SUMMARY 



1. Comparison 

By comparing Eq. (11.4) with the Blanchet and Faye 3PN equation of motion [40], we find the following relationship: 
miaf « work ^ rni{a^'^)^^_i^+mxnxU^,^+mra^\s^^^^, (12.1) 

where miO^'^'^ is the 3PN acceleration given in Eq. (11.4), (af ^)a=-i987/3080 is the Blanchet and Faye 3PN 
acceleration with A = —1987/3080, and miail^^j^^ is given in Eq. (11.2) with cRa replaced by for notational 
consistency with the Blanchet and Faye 3PN equation of motion shown in [40]. miaija^ is an acceleration induced 
by the following dipole moments of the stars: 

3709 

^i.BF = (12.2) 

We can compute Jnidij^^ Bp by substituting (5\bf instead of 5\q into Eq. (8.5). Thus, by choosing the dipole 
moments, 

D\&,BF = £*^A0 ~ £*^A,BF) (12.3) 

we have the 3PN equation of motion in completely the same form as (a?^)A=-i987/3080- In other words, our 3PN equa- 
tion of motion physically agrees with (af ^)a=-i987/3080 modulo the definition of the dipole moments (or equivalently, 
the coordinate transformation under the harmonic coordinate condition). In [47], we have shown some arguments 
that support this conclusion. 

The value of A that we found, A = —1987/3080, is perfectly consistent with the relation (1.1) and the result of [38] 

(^static = 0). 



2. Summary 

To deal with strongly self-gravitating objects such as neutron stars, we have used the surface integral approach with 
the strong field point-particle limit. The surface integral approach is achieved by using the local conservation of the 
energy momentum, which led us to the general form of the equation of motion which is expressed entirely in terms 
of surface integrals. The use of the strong field point-particle limit and the surface integral approach makes our 3PN 
equation of motion applicable to inspiraling compact binaries which consist of strongly self-gravitating regular stars 
(modulo the scalings imposed on the initial hypersurface) . Our 3PN equation of motion depends only on masses of 
the stars and is independent of their internal structure such as their density profiles or radii. Thus our result supports 
the strong equivalence principle up to 3PN order. 

The multipole moments previously defined in papers I and II are found to be unsatisfactory at 3PN order in the 
sense that these moments are not defined in a proper reference coordinate and consequently they contain monopole 
terms. By taking account of the effect of Lorentz contraction on these moments, we succeeded in extracting cleanly 
the monopole terms from these multipole moments up to the required order. 

At 3PN order, it does not seem possible to derive the field in a closed form. This is because not all the superpotentials 
required are available, and thus we could not evaluate all the Poisson-type N/B integrals. Some of the integrands 
allow us to derive superpotentials in series forms in the neighborhood of a star. For others, we have adopted an idea 
that Blanchet and Faye have used in [39-41]. The idea is that while abandoning complete derivation of the 3PN 
gravitational field throughout N/B, one exchanges the order of integration. We first evaluate the surface integrals 
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in the evolution equation for the energy of a star and the general form of the equation of motion, and then we 
evaluate the remaining volume integrals. Using these methods, we first derived the 3PN mass-energy relation and the 
momentum-velocity relation. The 3PN mass-energy relation admits a natural interpretation. We then evaluated the 
surface integrals in the general form of the equation of motion, and obtained an equation of motion up to 3PN order 
of accuracy. 

At 3PN order, our equation of motion contains logarithms of the body zone radii Ra- Practically, we cannot discard 
Ra dependences if Ra is in logarithms. We showed that we could remove the logarithms by suitable redefinition of 
the representative points of the stars. Thus we could transform our 3PN equation of motion into an unambiguous 
equation which does not contain any arbitrarily introduced free parameters. 

Our so-obtained 3PN equation of motion agrees physically (modulo a definition of the representative points of the 
stars) with the result derived by Blanchet and Faye [40] with A = —1987/3080, which is consistent with Eq. (1.1) 
and ecstatic = reported by Damour, Jaranowski, and Schafer [38]. This result indirectly supports the validity of the 
dimensional regularization in the ADM canonical approach in the ADMTT gauge. 

Blanchet and Faye [39,40] introduced four arbitrary parameters. In Hadamard's partie finie regularization, one 
has to introduce a sphere around each singular point (representing a point-mass) whose radius is a free parameter. 
In their framework, regularizations are employed in the evaluation of both a gravitational field having two singular 
points and two equations of motion. Since there is a priori no reason to expect that the spheres introduced for the 
evaluation of the field and the equations of motion coincide, there arise four arbitrary parameters. This is in contrast 
with our formalism where each body zone introduced in the evaluation of the field is inevitably the same as the body 
zone with which we defined the energy and the three-momentum of each star for which we derived our equation of 
motion. 

Actually, the redefinition of the representative points in our formalism corresponds to the gauge transformation 

in [40], and only two of the four parameters remain in [40]. Then they have used one of the remaining two free 
parameters to ensure the energy conservation and there remains only one arbitrary parameter A which they could not 
fix in their formalism. 

On the other hand, our 3PN equation of motion has no ambiguous parameter, admits conservation of an orbital 
energy of the binary system (when we neglect the 2.5 PN radiation reaction effect), and respects Lorentz invariance 
in the post-Newtonian perturbative sense. We emphasize that we do not need to a posteriori adjust some parameters 
to make our 3PN equation of motion satisfy the above three physical features. 

Finally, we note here that Blanchet et al. [43] have recently obtained the same value of A by computing a 3PN 
equation of motion in the harmonic gauge using the dimensional regularization. 
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APPENDIX A: Q^^^ AND E^^" 

In this section we evaluate Qa'^ and R^^^-' up to O(e^), the order required to compute the 3PN field. In the 

following, we shall omit the ei?^ dependence as explained in the appendix of paper II. In the evaluation of the field 
up to 3PN order, <io/i'^^ and <sh>^^, we use <s^^ as the integrand in the surface integrals of Q^'' and i?;^'*-* . Then 
straightforward calculation gives 

Qa' = e"' E ^" / ["^^" - nA]v1 Va'' = O(e') (Al) 

for ^ 1 and 
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„=4 J9B^ '- 



KA 
Vl 



_^4 3^^<ii> ^ Q^^5) for / = 0, 

£4 3^' J^^^fei^,^^ _ _ ^ij^k^^ ^ Q^^5) for ; = 1^ (A2) 

O(e^) otherwise. 

ii^'*^' contribute to the field at 3PN order through the moments Z^'''^ . See Eqs. (3.21)-(3.30). Combining the 
above results, the Q^'^ and i?^'*-' contributions to the 3PN field /iq^ become 



A=l,2 ''^ ^ / A=l,2 

with Cqa = 12 and /ig^ = O(e^), where we used Eqs. (A8), (A9), and (AlO) below. 
The hq'i^ field affects the equation of motion, 

'12 '12 

where a\Qj^ is the QR field contribution to the acceleration of the star 1. Note that in the above equation we have 

not yet taken into account the effect of the 3PN Q\ integral {eQ\)- gQa does not affect the 3PN field, but the 3PN 
equation of motion through the 3PN momentum-velocity relation. We evaluate eQ^ VIII and Appendix E. 

For convenience, we list Q^'* integrals, R^'^-' integrals, and the body zone contribution including multipole moments 
up to 3PN order. We retain here (and only here) the dipole moments of order O(e^), since it is appropriate to define 
the representative points of the stars by D\ = e'^M^v\ + O(e^) when we are concerned with the spin effects (see 



2m2Mf 



paper I). 

Qi=^' ""^Z2 •"^ +0(e^), (A5) 

and Qf'' = 0(e5) for / 7^ 0, 



of 



12 



g^3^ (3/i^' - Wnl'nl^ni^ - 12/f nfanla + 3/^71^2^2 + 3/^ni2ni2) 

42^2 / , k[lk\l k J r,k[H]l k 3 y^ml 

'^^ \ l^Pl'i^l ~'^12"'12^1 ^ "-12"-12^1 

- Zf'^"" - + (monopole part) + 0{e'), (A6) 

^ ^4^ (2<5*0/^)'nl2 - 2/i(^nt^ - PHi^^^^ + 2,Ii^^^^ 



12 

- ,5*'=zf"'l"n'i2 + Zf'^'n\^ + zf^^'v?^^ + "l'nf2 + zf'^^n?^^ 

- 3Zf["=l'ni2 - 3Zf ''^^ nl2) + (monopole part) + 0(i), (A7) 



and Ry''^ = 0{e^) for Z 7^ 0, 1. 
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In the following, we have used Eqs. (3.19)-(3.30): 
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(AlO) 



We restrict our attention to the equation of motion for two spherical compact stars in this paper, however our 
formulation can be extended to an extended body with higher multipole moments, as shown in paper I. 



APPENDIX B: DERIVATION OF EQ. (3.31) 

In this section, we show a derivation of Eq. (3.31) without using the Dirac delta distribution. The following proof 
is essentially due to [60]. 

Suppose that a two-dimensional surface S surrounds a three-dimensional volume V. Suppose that there is a point 
P{x) and a point Q{y) both inside of V, and the distance between the two points is r. Thus, 

r = \x — y\. 

Note that \7'^{l/r) = except for r = 0. 

Define a sphere V which is centered at P{x) and has a radius d that is sufficiently small so that V is enclosed 

completely by V. The surface S of y divides V into two regions. We call the outer region V". 
The Green's theorem (e.g., [61]) states for a certain function v — v{x) that 

Note that V" docs not include the point P{x). Wc assume here that for v{x) the second derivative with respect to x 
exists and is continuous in V" . In the end of this section, we mention whether this assumption holds for our particular 
example, Eq. (3.31). 
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The third term can be evaluated as 



where rij is the outward normal of the surface E. If dv{y)/dy^ is finite in the d ^ limit, then the first term is zero 
in this limit. Thus, 

'''^'^ viy)^('-)]dS. = 4Mx^), (B3) 



in the — > limit. 
In the same manner, 



Then we have 



r dy'^ dy^ 



^-^d^y=liinj V^viy)rdrdn = 0. (B4) 



d^y+f\ -^rr -<y)in\-]\ dSi + 4.-kv{x) = o. (B5) 



r dy^ dy 

This is essentially Eq. (3.31). 

In the particular case of Eq. (3.31), V = N/B. We recall that the stars are completely enclosed by the body 
zones Ba and thus there exists no matter in N/B. The integrand v{x) or f{x) and g(x) in Eq. (3.31) have thus no 
singularity in N/B even in the point-particle limit. These functions have a smooth second derivative with respect to 
X in N/B. 



APPENDIX C: CORRECTION TO THE MULTIPOLE MOMENTS 



A natural reference coordinate where we would define multipole moments of a star may be a coordinate in which 

effects of its orbital motion and the companion star arc removed (modulo, namely, the tidal effect). In other words, 
such a natural reference coordinate may be the generalized Fermi coordinate [59], in which the metric is Lorentzian 
at z\{t). We define the two stars to be spherical in such a coordinate in this paper [30]. 

In paper II, wc defined a multipole moment, which we call the NZC moment in this section, as a volume integral 
over the body zone which is spherical in the near zone coordinate (NZC). Then a question specific to our formalism 
is whether the NZC moments affect the orbital motion of the stars which are spherical in the generalized Fermi 
coordinate (GFC). 

At lowest order (0(e°)), the NZC multipole moments and the GFC multipole moments defined as volume integrals 
over a sphere in GFC must be the same. For a spherically symmetric star, the NZC moments and the GFC moments 
with trace-free operation or antisymmetrization on their indexes vanish at the lowest order. We further assume that 
the GFC moments and the NZC moments are zero for a spherical compact star at the lowest order because of its 

compactness. 

For the mass monopole, we have already obtained the relation between the NZC monopole (the energy) and the 
mass via the evolution equation of the energy. Thus, we seek the IPN correction to the NZC Lth multipole moments 
with L > 2 since physically relevant multipole moments appear at 2PN order and we are concerned with the 3PN 
equation of motion in this paper. In the following, we compute the corrections which do not include higher order 
multipole moments than the energy monopole. 

As for construction of GFC, we assume, up to the relevant order here, that the transformation from NZC to GFC 
or vice versa takes apparently the same form for a strongly self-gravitating star in which we are interested here and 
for a weakly self-gravitating star for which GFC has been constructed in [59] . This is because the construction of the 
generalized Fermi coordinate depends on how the star moves, and because the equations of motion for binary stars 
take the same form regardless of the strength of the stars' internal gravity up to 2.5PN order (paper II). An important 
difference is that the mass parameter in the transformation from NZC to GFC for the strongly self-gravitating star 
includes the strong field effect, as explained below Eq. (3.36). 

We now assume the following coordinate transformations from GFC (f, rr'^) to NZC (r, a;*'): 

y' = z\{T)+y\ + Sy\T,y'X), (CI) 
r-rp = f-Tp + 6T{T,y\), (C2) 
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with 



n=l 



(C3) 

(C4) 



ra=0 



where the capital index denotes a set of collective indexes, e.g., A*"^" = A'-J'i = — z\. In Eq. (CI), y\ is on 
the f = const surface, rp and fp are fiducial time coordinates. We assume that A^-^^jB"^" = 0{e^). 

The NZC moments are defined on Ba{tp) which is a sphere centered at {tp,z\{tp)) with radius eRa in NZC, while 
the GFC moments are defined on Ba{jp) which is centered at the same event {tp,Za) in GFC {{tp,z\{tp)) in NZC) 
with the radius cRa in GFC. The corrections we shall compute are then 



,2;+4— s jJil^i' 



^ -'A,GFC' 



JBa(,t=tp) 

d'yAf^A^J,{r,y% 



21+4,— s jJittv 

-"a.gfc 



(C5) 

(C6) 

(C7) 



i(f=f_p) 



where s = 2 for (/i, t^) = {hi) or otherwise. K'^'',{f.ff) = A'^''{T,y\). 

We now express Ianzc using the generalized Fermi coordinate. Note that y^'(T) in Eq. (C6) is on the r — rp — 
const surface. A necessary coordinate transformation that relates y\(T) on the r = const surface with y\ on the f = 
const surface can be obtained by modifying Eq. (CI) using retardation expansion. Up to IPN order, the result is 



y' = z\{T)+y\-5y\Tp,y'X), 



with 



Sy\rp,y'X) = ^ {B-'-irpyAirp) - A'^"(tp)) y^". 

n=l 

Using tetrad e''a(f, y') = 5''^ + O(e^), we have 



(C8) 
(C9) 



,2i+4-s jJmv 

E -'a.nzc 



BA{f=f(fp,y\)) 



d-'yA 



d{y\) 



I 



y^Y{[y'X-5y^'{r,y'A)} 

fc=i 

xe 



^^{T,y\)e-^{TSA)K-{r,f) 
d'yAyjlA'S',{r,n 



(CIO) 



where in the last equality we neglected all the corrections that result in the multipole moments of order L' > L. The 
integral region Ba {f — T{Tp.y\)) which corresponds to Ba{t = tp = const) is not a f = fp const 3-surface nor 
spherical in GFC. In fact, from Eq. (C2), 



Ht^p^Va) = rp- 5T{Tp,y\)- 



(Gil) 



Thus we make a rcitardation expansion around f = fp in the integrand. Then to make the slightly perturbed sphere 
Ba (f = fp) into the sphere Ba (f = fp) with radius cRa, we change the integration variable y\ into y\. Up to IPN 
order, the transformation may be obtained by inverting Eq. (C8), 



f = ZA{Tp) + fA + Sv\Tp,rA)■ 
\^S\ng Eq. (Gil) and Eq. (C12), we simplify Eq. (CIO) up to IPN order as 



(C12) 



42 



,2;+4-s jJilJ-v 

^ -'a.nzc 



BA{fp) 
JBAifp) 



BA{fp) 



d{y'A) 



d 



x{K'^^^ [rp,y\ + Sy\Tp, y^) ) " Srirp, f^)—A^J,{rp, f) 



Ba{tp) 



y^^6T{rp,y\)—A^J,{fp,f)\. 



(C13) 



As before, we have discarded terms which end up with multipole moments of order L' > L. Integrating by parts and 
rewriting by y^, we finally obtain a formula for 51^^'^ , 

= J2 (B^"(Tp)^;3'(rp) - A-^"(rp)) / d5„yi'^f"A^';(fp,r) 

1 JdBA(fp) 



d 



n=l 



IdBAifp) 



P JBAifp) 



d'yAyiyr^^J'irp,f) 



5^(B^"(rp>3'(rp)-A™^"(rp)) 



n=l 



BA{fp) 



d'yA^yiA'J'{rp,n 



(C14) 



Notice that dz\{Tp) / dip = 0. 

The last two volume integrals in the above equation result in the multipole moments and we neglect them here. As 

for the first term, we can evaluate the surface integral to IPN order explicitly by replacing all the hatted quantities 
by those without a hat (and G' by N). We foimd that ^A^'' (n < 7) do not contribute for any L. On the other hand, 
8 Ay for L = 2 and n = 1 in the summation in the first term of Eq. (C14) gives a monopole correction to the 3PN 
gravitational field (through the IPN correction to the quadrupole moment which itself appears at 2PN order in the 
field), 



dBAifp) 



dSmyAyAyAAG'{rp,y") = 



dBA{rp) 



dSMAyAA]^ijpMl + 0{e')) 



dS^y\y'AyAsANirp,f)+0{e'^) 



10\ 



dBAirp) 

,4m^ 



Neither s^'^ for L > 3 nor n > 2 contribute to the 3PN field. 

The coeflacients A'"'=(rp) and -B'=(tp) may be read off from the results in [59,62], which are 



B\Tp)vi{Tp) = e\\vi+0{€ 

,2 



A^i(rp) = e'[lvlvi-'^6^^]+0{e') 
'2 ri2 



(C15) 



(C16) 
(C17) 



for the star A = 1. Thus, the coefficient in front of the surface integral, Eq. (C15), becomes e^{v\vl/2 + 6^^m2/ri2)- 
Finally using the coefficients above, we obtain the IPN corrections of the multipole moments that affect a 3PN 
equation of motion for spherical stars (we rename (57^^^ by SP^), 



43 



Sr^ = ^ {-^^v\v\ - - ^^^A + 0{e% (C18) 

\ 5 5 ri2 J 

A similar equation holds for the star A = 2 by exchanging 1 <-+ 2 in the above equation. Notice that since the 
quadrupole moments at 2PN order is symmetric-trace-free, the last two terms do not contribute to the 3PN field. 

The correction SI^^-'^ should appear even when m-z 0. Not surprisingly, with the correction Slf^-'^ , the 3PN 
gravitational field for a single star moving at a constant velocity derived by solving the harmonically relaxed Einstein 
equations iteratively agrees with the boosted Schwarzschild metric in the harmonic coordinate up to 3PN order. 

It is possible to use Eq. (C14) for the L =1 case and the result becomes again the 3PN correction to the field (we 
again rename 51'^'^ as 6D\), 

6D{ = e^^ri,. (C19) 

(In the computation of SD\, the surface integral with sXn"^ ,a0 as the integrand is found to vanish.) However, any 
change of the dipole moment amounts merely to a redefinition of the representative point of the star and it causes no 
physical effect. In fact, we have not taken into account SD\ to derive our 3PN equation of motion, Eq. (11.4). 

APPENDIX D: RENORMALIZATION OF THE MULTIPOLE MOMENTS 

This section explains the "renormalization" of the multipole moments and that the field does not depend on cRa- 
This is rather trivial, however we show this section for clarity. 

We first define a symbol part^^j^ F which is the ei?^-dependent part in an expression F except for the logarithmic 
dependence of cRa- Correspondingly, we define disc^fl^ F, which means to discard all the eiJ^-dependent terms in F 
other than Inei?^. By construction, 

F = discF + p&vtF. (Dl) 

We give an example of disCeK^ in Sec. VC. 

Now, to derive the field, we study the following Poisson integral for a certain function f{x). f{x) is some combination 
of A''^(t, .f) and is assumed to be nonsingular in N. For notational simplicity, we do not write time dependence 
explicitly in f{x), 

Jn f - y\ ^2 ■'ba f - yl Jn/b \x - y\ 

The second volume integral is evaluated with the help of the superpotential g{x) that satisfies Ag{x) = f{y) in 
N/B. Using Eq. (3.31) and expanding the kernel l/\x — y\ around yA = V ^ ZA^ we obtain 

[ -fy- fU-f^ - \- \- (2n-l)!!r<^"> /' j,^ dg{yA + ZA) 



+ E E ^ Jn.i i dS,g{yA + ZA)f 
+ •••. (D3) 

Here we only show explicitly the terms which possibly depend on cRa, and "• • • " denotes ei^^-independent terms. 
Thus, using the symbols introduced above, we have 

/ rS^m=dmc[ -^/(^+part/ jS^fiy), (D4) 
Jn/b \x-y\ <^Ra Jm/b \x-y\ cRa Jn/b \x - y\ 

with 
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part 



cRa Jn/b \x-y\ 



' part 



EE 



(2n- l)!!r^ 



A=l,2n=0 

(2n- l)!!r^" 

A=l,2n=0 



— 

Ba "Va 



za) 



+ EE' 



n\r 



dSkgivA + ZA 



Ba 



(D5) 



For the first volume integral in Eq. (D2), we use the multipole expansion 



E / w^f^^ 



J2 ^(2n-l)!!r<^">/^ 



A=l,2 n=0 



with 



7^"= / d'yAf{yA + ZA)y^''. 

Jba 



(D6) 



(D7) 



7^" is the multipole moment of the star A. We simplify here the definition of the multipole moments so that we 
omit the scalings on the integrand A^'^(r, x) and hence f{x). Obviously, this definition is enough to study the cRa 
dependence in the multipole moments and the field. 

7^" in general depends on cRa since the integrand f{y) is noncompact support. A possible cRa dependence in 
I^" may be examined by the following. First, since we are studying the nonsingular sources, we expect that f{x) is 
smooth so that we can assume Ag{x) = f{x) just inside of OBa- Thus, the cRa dependence in 7^" can be examined 
via 



rKr, 



in the neighborhood of dB^ 

dSk 



dBA 



K„ dgiVA + za) 
Va 



dVA 



d^yAAg{yA + ZA)yA" + ' 

^ givA + za) 



dy^- 
dyA 



Here "• • • " again denotes terms that do not depend on cRa- 
Using the symbols introduced in this section, we have 



7^" = disc 7 

eRA 



part 71"", 

eRA 



disc 7"" 

eRA 

part7j" 

eRA 



a"' 



disc / d^yAfivA + ZA)y 

JBa 

part / d^yAf{yA + ZA)yA 

eRA J Ba 



part 

eRA 



,o ( K„dg{yA + ZA) dy^ ,^ ^ ^ . 

Va ^nt jnrdiyA + ZA) 

dBA \ dy'X dyl 



(D8) 

(D9) 
(DIO) 

(Dll) 



We call 7^^ = disceij^ 7^" the "renormalized" multipole moments that are independent of powers of cRa by definition 

of disc^fl^.' 

It is clear that the cRa dependences in the multipole moments Eq. (Dll) cancel out those of the N/B contribution 
shown as the first and the second terms in Eq. (D3). In conclusion, we find that we can discard (or neglect) all the 
e7?^-independent terms other than In cRa terms when we compute the field, 
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Jjv F - y] ^2 ■'Ba f - v\ Jn/b \x-y\ 



IN/B 

(2n-l)!!r<^">/J" ^ f d?yf{y) 



^ ^ (2n-l)!!rr""/r . / 



> > T-ZT part IT + part / , j^"^ 



In the main body of this paper other than this section, we write as and omit the symbol disc^ij^ in front of 
the Poisson integral over N / B for notational simplieity and from triviality of the fact that the total field is independent 
of ei?^. As an exception, we write disCe/{^ in Sees. V and IX to make our discarding cRa procedure clear. 

Finally, we mention that it is straightforward to extend the arguments here to show that the cancellation of cRa 
terms between the body zone contribution and the N/B contribution occurs for any retarded field, that is, n > 1 
terms in Eq. (3.10). 

APPENDIX E: x PART 

We derive the functional expressions of P^^ on tua, v\, and r^2- Here we defined as 



[ d%xr^,ap. (El) 

JBa 



By the definition of Xjv'"'^>a/3) 



Ba 



kl, 



kl, 



thus we can obtain the functional expressions of P^^ using Gauss' law. In fact, up to 3PN order, the definition of 

Pax gi^^s 



4y2 ^ rn2 2m, 



ri2 ri2 



ri2 



2mf 5m-]mo mi mi / 14 „ 11 o 



22 2 20 

- y (^1 • ^2) - 5(^12 • ^1)^ + y ("12 ■ 'i^l)(^12 • vi) - 'i{ni2 ■ 

ma / 197 2 38 _ ^, 2_ ^ 2_ ^, 

+ + —V2 - —{Vi ■ V2) + —("12 • Vi) - -{ni2 ■ Vi){ni2 ■ V2) 

ri2 \ 60 6 6 15 3 

224 34 2 2 44 44 2,-, -,^ 82,-, ^^ 8_ ^,2 

+ Y^^^l + Y5^1^2 + 3^^2 - Y^MVI ■ V2) - -^V^ivi ■ V2) + —{Vl ■ V2) 

- ^^^1(^2 • V2f - ^vl{ni2 ■ V2)^ + ^('^1 ■ '^2)(ni2 • V2f 

+ 0{e'). (E2) 
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PJ^ of O(e^) affects a 3PN equation of motion only through the field loh'^'^. 

The dipole moment and Q\ integral of the x P&^t give a nonzero contribution starting from 3PN order in the 
momentum-velocity relation, 

Dl, = .''-I^ri,+Oi.% (E3) 

«i« = (^JJnSlA + „„f + oi.^y (E4) 

Here again, we used Gauss law to derive D\^. On the other hand, we evaluate directly from the definition of the 
using Gauss law and found that 

PX^ = PX^v\ + Q\^ + e'^ + Oie') (E5) 

is an identity up to 3PN order. Thus, the representative point of the star is defined with the part of the momentum- 
velocity relation, not with the x part. 

Finally, by evaluating the surface integrals in the evolution equation, 

^ = -e-'<f dSuxT^,a0 + e-'v\i dSkxT^,a^, (E6) 

JdBA JoBa 

we found that the resulting equations for dP^^/dr are consistent with the explicit expressions of P^^ directly obtained 
from their definitions, Eqs. (E2) and (E5), as expected. 



APPENDIX F: LANDAU-LIFSHITZ PSEUDOTENSOR EXPANDED IN EPSILON 

The Landau-Lifshitz pseudotensor [53] in terms of h^"^ which satisfies the harmonic condition is as follows. 

+ \ {g'^'^a"^ - Ig'^-'g"^^ (g-^sgec - lg-yeg6?j h?\<.h'^,p- (fi) 

We expand the deviation field h'^" in a power series of e; 



n=0 

Using this equation, we expand in e. 

Here wc only show io[— IGTrgt^"^]. See paper II for <9[— IGTriji^'^]. Note that all the divergence such as h^^ in 
paper II should be replaced by —Ji'^'^^t consistent with the following results. This is simply because it is practically 
much easier to use h^'^^r than —h^^'^^k, 

,o[-lQngtll] 

7 1 

= — ^4/i^^,fe8/i^^' + ^4^^^' eh k,i — 4/1^^, fee^'" ,T + 24/1'^ ' qK^ (k,i) 

3 7 1 

7 1 

— g6^^^,fe6/i^^' + ^ih'^^ ,Tih k,T + ih'^ ' 4hki,T 

-|--4'l Ankl^m — ■^411 4n'km,l — g4'l fe,m4'l I 

17 7 7 

+ -4^^ 4h^^ ,T4h'^^ ,k + ^4h 4/1^^,^4/1^^,; — g (4/1^^) 4h^^,k4h^^' + ^6h^^4h^^ ,k4h^^' 

-24h^>'ih^\uh^\k - |4/i"\/i""''4/i"fc,i, (F2) 
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io[-167r5ili] 

4 ' 4 ' 

-^4/i^^4/i^^r6/i^^'' + \Ah^\h^^'^QK'\k - \Ah^ kih'^'^Uh^^^^'^^ + 24h^^ ^k^^^'^^h^' 

-^4/l"(4/l"",r)' - ^4/l''=4/i"",r4/l"",fc + ^ (4/l"")'4/l"",r4/l""'' - ^6/l""4/i"",r4/i""'' 

4 z 

— -^iah'^^Ah'^^,kAh'^'^'^ + 24^^'°4^^fc,T4/l^^'* — ^h^^ dl^'^ ,uh'^\k + 2(4/l'^'^)^4/l'^'^^fe4ft,'^['^'''' 

— 26/l'^'^4/l^'",fe4/l'"''^''' + 4/l"4/l'"*''4/l'"[;,fe] + 24'^^'s4/l^^,r4/l'"'^'' 

+24/l^\/l^fc,i4/l"['''l +24/l"S/l"/,M/i"l'''> +4/l"\/i"^r4/l"\fc +4/l'='4/l"",fc4/l",i 
+ j4/l"4/l""''=4/l';,fe - ^4/i"fe4/lV''=4/ll""l''^ + 24/l"fe4/i"",M/l'='''^ 
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io[-167rfft^J 

+ \{S\S^i + ShS'i - 6'' Ski) ieh^^^" + iK^m^'' - 24/1^^4/1""'' + ^Ah^\r) sh^^'' 
+liS'kS^i + PkS'i - S'^Ski) ieh^^''' - ihJ^m'^ - 4/i"%/i""'*=) ^hTj 

+24/l'=[''\/l^fc.i+24/l'=l^'\/i'fc,i 

+ {S'kPl + 5h6S - S'^Skl) (^ih''\h7m'^Qh7"''^ - 4/l""4/^"^m6/l"'"'' + ^4/l"m4/l""''=6 
+24/l""4/l"['''=l6/l"V + 24/l""4/l"[^''=l6/l",/c 

+ ^S'^ (24/l«4/i^'=''"4/l^',m - 24h''Uh^"'Mh^m,i+ih''Uh^\k6h^\i - Ah'^h'^Uh^^kAh^^ 

+ l4h'^Uh^\r f - ^4h^^4h^hh^\k4h^^''' + 4h^\h^\k4h^^^^Uh^^'' +4h'hh''\k4h^\r 

o 4 

—4h^^4h'^'''''4h^[kjj + 4h^\h'^\k4h'^'' ,1 — 4hki4h^*''^ 4h'^^'^ 



TT 

k 



+ \4h'' 4h^\kf^h7^-'' - \4h:'\k&h7^-^\h^^'' - }-^h-\^^h---(\hi)k 

+5^^ ^24/l'^fe4/l'^;_TO4/l*[''"'l — -4h'^ k4h^'''Uh"^ m,l — -^{4h'^'^)'^4h^^''' 4h^ l,k 
+ ^4/l^'^4^^^,r4/l'i,fe 
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+4/l"fe4/lV^'4/il"'=l'^'' + 24/l"'=4/l^['''l4/l«'^' + 24/l^*=4/l"[^'']4/lfe;'' 
+24h^>'4h^^^''hhh,l + 24/l"\/l"I''^'l4/i'fc,; - ^4/i""4/l'=fe'^'6/ll""l'^'' 

3 11 1 

— -4/l'^'^6'l^^4'l^^,fc4^^^' + ^8'l^^4ft-^^,fc4ft-^^' + -^ih'^ kih'^ 4h'^'^ ,l4h^^' — 2^^^ 4'i^^,t4^^/c,t 

— [Ah'^'^Y Ah'^^''^ ih'^ k,T + &h'^'^ Ah'^'^'^ Ah'^ k,T + Ah'^'^ Ah'^^ Ah^'^'\h'^ k,l — Ah'^ kAh'^^ ,lAh'^^ ,T 
+ {ih'^^Y ihJ^'^h'^ [k,l\ — 6h'^^Ah^'''''4h'^[k,l] — 4/l^'^4/l'"z,fe4/l'"',r + ^4/l*^Z,T4ft'/c,T 

+ -^4h^^ 4h'^^ ,k4^"^ m,l — -^4h^k,T4h\^T + -^4h^ k,T&h'^^ ,t — -^4h''^ 4h'^^'"^4hkl,m + ^4/l'"^4/l'"^,r6^'"^,T 

1 3 1 

+ -Ah^ 4/l^^,fe6/l^^,T — o(6/i^^,t) +74^^" 4^^^,r6'l^^,fc + 4^^^4/l'" ,TQh'^'^ ,k 

4 o 4 

— -^Ah^^ 4h^ k,iQh^'^'^ — -(4/1^^)^4/1^^, fce/*'^^''' + 2^}^^^ Ah^^ ,k6h^^''^ + -^Ah^^eh'^^ ,k6h^^'''^ 

+4/l^''4/l^'',r4/l^^'^'4/l^'^^ - (4/l^'')^4/l^^''4/l^^"'' - ^4/l^''4/l'"fe4/i^^''4/l^^"'' 

+ ^4/;,^^6/;,^^4/l^^'*4/l^^'j' - i8/l^^4/l^^'*4/l^^'^' 
-24/l^^4/i^fe4/l^^'^*4/l'^''''^'^ + 24/l^fe,r4/l^'''^*4/l^'^^ 

+2(4/l^^)24/l^fc'^4/i^f''''l - 2(4/l^^)^/l^^fe4/l^f''''l + 2(4/l^^)\/l^^'^'4/l^'^\r 
-26/l^^4/l^fc'^4ft^[''''l + 26/1^^4/1^-'', fc4/l^'''''l - 2eh''\h''^'^\y'^-' ,r 

-2Ah-'\rAh^^'Ah^^^ ,r + ^aK" kAh'''''^\h^^\r - 2Ah^'' ,rAh^^' aH^^^ ,k + 2Ah^ kAh^^''^''^ Ah^^^ ,r 

+ \Ah}' Ah^^-\h\^k - \Ah'ukAh^^^^\y^'' - ^Ah^\kAh'r^\h^^'' + ^Ahkuh^^'^Uh^"'^'^^ 

-4h'\r4hh,r - \6h''\r4h^^'^\h'^^ - ^4h^\r6h^^'^Uh^^^ + ^{Ah^^fAh^^''-'eh\^^^'^'> 

-6/l^^6/l'"'"'^'4/ll'"'"l"''^ + 4/l^fe4/l'"^'^'6/l'^^''^'' - 24/l^^6/l^^'^*4/l''''^,r - ^4/l^'"6/l'"'"''6/l^'"'^'- (F4) 



APPENDIX G: SUPERPOTENTIALS 

Here we list useful superpotentials. Below, satisfies A/^™'") = rj"r^. Similarly, A/f™-'"'"''") = 

ri"r2 Inri lnr2. The appendix in [33] gives a greatly useful list of superpotentials including /^''^'^\ f^^'^\ /'^~^'^\ and 
j"(3>-i). Other useful superpotentials are given in [28,40], 
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^(-2,2|^!l^l,_^j+2r;,+r?)l„n, 



rfan \rij r^2f2 \r2 J rf2rir2 



rm2 \r\ 



& 



3,0) 



Inn 

r-1 ^ 



/(-3,i) = _!i + i„5 + !iai^fl 



/ 



(-4,-1) 



?'2 



2^2 2 ' 
^'l'l2 



y(0;ln,0)^ 1 5 2_ 

6 ^ 36 ^ 

An example below shows how we construct superpotentials from simpler ones, 

j(6,-2) _ _g^(0;0,4:ln) _ 3^(2,2) _j_ ^^r^a y.(0;0,2an) _ '^'^^12 j(2;0,0;ln) _|_ 6^12 j(4,-2) 
1^1 '^24i*^24 ^2i^li *^42i ,'^24i 

14 " 35''^^''^ TO''^^''^ ~ 49 y ^''^ ~ ''2 111 ''2 + ^rirj Inra, 
which is computed from /(0;0,4;i„)^ j(2,2) j(o;o,2;in) j(2;0,0;in)^ ^^^^ "2). They are 

.(0;0,4;ln) _ J_„6 1 _ „6 

/•(2,2) _ _!i , J_ 4 2 , J_ 4 2 



^(0;0,2;ln)^ 1 4i^^^_ 9 



20'2^ 400'2' 

y(2;0,0;in) ^ _ JL^4 ^ J_^2^2 _ JL^2 2 ^ J_^4 ^ Ir^r^lnrz + — r^r^ In r2 - — r^lnrs, 
•' 200 1 60 1 12 ;^go 12 2 -r gQ 2 T -j^Q 1 2 2 T -j^g 12 2 2 20 2 2, 
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